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ABSTRACT 

We present results of numerical simulations which examine the dynamical stability of known planetary systems, 
a star with two or more planets. First we vary the initial conditions of each system based on observational data. 
We then determine regions of phase space which produce stable planetary configurations. For each system we 
perform 1000 ~ 10 6 year integrations. We examine v And, HD83443, GJ876, HD82943, 47UMa, HD168443, and 
the solar system. We find that the resonant systems, 2 planets in a first order mean motion resonance, (HD82943 
and GJ876) have very narrow zones of stability. The interacting systems, not in first order resonance, but able 
to perturb each other (v And, 47UMa, and SS) have broad regions of stability. The separated systems, 2 planets 
beyond 10:1 resonance, (we only examine HD83443 and HD 168443) are fully stable. Furthermore we find that the 
best fits to the interacting and resonant systems place them very close to unstable regions. The boundary in phase 
space between stability and instability depends strongly on the eccentricities, and (if applicable) the proximity of 
the system to perfect resonance. In addition to 10 6 year integrations, we also examined stability on ~ 10 8 year 
timescales. For each system we ran ~ 10 long term simulations, and find that the Keplerian fits to these systems 
all contain configurations which may be regular on this timescale. 
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1. INTRODUCTION 

By September 2003, 13 planetary systems had been discov- 
ered (including our own Solar System). The extra-solar plane- 
tary systems (ESPS) are, possibly due to observational biases, 
markedly different from our own in several ways. In our solar 
system (SS) the Jovian mass planets all orbit at distances larger 
than 5 AU, and on nearly circular orbits (e < 0.05). ESPS, on the 
other hand, contain giant planets in a wide range of distances 
and eccentricities; some are 10 times closer to their primary 
than Mercury, and others orbit with eccentricities larger than 
0.5. In this paper we attempt to categorize these systems dy- 
namically, constrain the errors of the orbital parameters, com- 
pare our SS to ESPS, explore the long term stability of each 
planetary system, and determine the mechanism(s) that main- 
tain stability. 

We examine these systems through numerical simulations. 
The integrations begin with slightly different initial conditions 
in order to probe observationally allowed configurations. This 
exploration of parameter space permits a quantitative measure 
of the stability of each system, and, hence, predicts which dis- 
tribution of orbital elements will most likely result in a stable 
system. In addition, a comparison of stability between systems 
may reveal which elements are most critical to the stability of 
planetary systems in general. 

A quick inspection of the known systems reveals three ob- 
vious morphological classifications: resonant, interacting, and 
separated. Resonant systems contain two planets that occupy 
orbits very close to a 2:1 mean motion resonance. A 3:1 sys- 
tem, 55Cnc (Marcy et al. 2002), has been announced, but its 
dynamics will not be examined here. Interacting systems con- 
tain planets which are not in mean motion resonance, but are 
separated by less than a 10:1 ratio in orbital period. These sys- 
tems are not dynamically locked, but the planets perturb each 
other. The SS falls into this category. The final classification is 



systems in which the (detected) planets' orbits are beyond the 
10:1 resonance. These planets are most likely dynamically de- 
coupled, however some of these systems warrant investigation. 

We examine planetary systems on two different timescales. 
First, we explore parameter space in 10 6 year integrations. For 
these simulations we vary initial conditions to determine stable 
regions within the observed errors. Second, we continue several 
stable simulations for an additional 10 8 -10 9 years. From these 
runs we then learn how robust the predicted stable regions are, 
and we also determine the mechanisms that lead to stability. 
Specifically we evaluate the hypothesis that some stable system 
require secular resonance locking. 

In many ways this paper performs the direct analysis that is 
approximated by MEGNO (Cincotta & Simo 2000, see also 
Robutel & Laskar 2001, Michtchenko & Ferraz-Mello 2001, 
and Gozdziewski 2002). MEGNO searches parameter space 
for chaotic and periodic regions. Our simulations show that 
to first order MEGNO's results do uncover unstable regions. 
In general, however, chaotic systems can be stable for at least 
10 9 years, as shown below. Our SS also shows chaos on all 
timescales (for a complete review, see Lecar, et al. 2001), there- 
fore only direct N-body integration can determine stability. 

This work represents the largest coherent study of planetary 
system dynamics to date. Our simulations show that the true 
configurations of most planetary systems are constrained by 
just a few orbital elements (or ratios of elements), and that sta- 
ble regions can be identified with integrations on the order of 
10 6 years. We also find that stability, as well as constraints on 
stability, are correlated with morphology. Resonant system sta- 
bility depends strongly on the ratio of the periods, interacting 
systems on eccentricities, and separated systems are stable. 

This paper is structured as follows. In §2 we describe the 
generation of initial conditions, integration technique, and in- 
troduce the concept of a stability map. In §3 we analyze the 
results of the resonant systems HD82943 and GJ876. In §4 
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we examine the interacting systems of 47UMa, v And, and the 
SS. In §5 the separated systems, specifically HD168443 and 
HD83443, are discussed briefly. In §6 we summarize the re- 
sults of §3-5, In §7 we discuss possible formation scenarios and 
inconsistencies between this work and other work on planetary 
systems. Finally in §8 we draw general conclusions, and sug- 
gest directions for future work. 

2. NUMERICAL METHODS 

In this section we outline the techniques used to perform and 
analyze these simulations. This paper follows the example of 
Barnes and Quinn (2001, hereafter Paper I). First we describe 
how the initial conditions of each of the short term simulations 
are determined. In §2.2 we describe our integration method. 
Finally in §2.3 we describe a convenient way to visualize the 
results of these simulations, the stability map. 

2.1. Initial Conditions 

In order to explore all of parameter space we must vary 5 or- 
bital elements, (the period P, the eccentricity e, the longitude 
of ascending node f2, the longitude of periastron, 03, and the 
inclination i) and every mass in the system. The period and the 
semi-major axis a are related by Kepler's third law. The argu- 
ment of pericenter, ui is the difference 03 and 0. 

For each system we perform 1000 simulations, each with dif- 
ferent initial conditions. Orbital elements that are easily deter- 
mined via the Doppler method (P,e,03) are varied about a Gaus- 
sian centered on the nominal value, with a standard deviation 
equal to the published error. We do not permit any element to 
be more than 5a from the mean. For the i and il elements, flat 
distributions in the ranges 0° < i < 5° and < £1 < 2tt were 
used. Note that this distribution of i and permits a maximum 
mutual inclination of 10°. These randomized orbital elements 
are relative to the fundamental plane. 

The Doppler method of detection does not produce a normal 
error distribution. As Konacki & Maciejewski (1999) show, the 
error in eccentricity can have a large tail toward unity. However 
their method, or other statistical methods, such as bootstrap- 
ping, require all the observational data (including reflex veloc- 
ity errors) to determine the shape of this error curve. As not all 
the observations of multiplanet systems have been published we 
are forced to use a normal distribution in order for comparison 
between systems to be meaningful. We therefore encourage 
all the observational data to be published, as some of the re- 
sults presented here (specifically comparing the percentage of 
simulations that survived) are less meaningful because we are 
unable to accurately estimate the error distributions of the or- 
bital elements. For completeness we also vary the primaries' 
masses, M„, based on other observations (generally determined 
via spectral fitting). However as the stars are all at least 100 
times more massive than their planetary companions, the slight 
variations in primary mass should not affect stability. 

The masses of the planets are determined by the following 
relation 



\sin(cos l (sini)cosQ)\ 

where m is the true mass of the planet, and m /,, is the observed 
minimum mass. By varying the mass this way, we account for 
all possible orientations, and connect the inclination of the sys- 
tem in its fundamental plane, to its inclination in the sky. Note 

'SWIFT is publicly available at http://k2.space.swri.edu/~hal/swift.html 



that this scheme requires the azimuthal angles in the planetary 
systems to be measured from the line of sight. 

2.2. Integration 

With initial conditions determined the systems were then 
integrated with the RMVS3 code from the SWIFT suite of 
programs 1 (Levison & Duncan 1994). This code uses a sym- 
plectic integration scheme to minimize long term errors, as 
well as regularization to handle close approaches. The initial 
timestep, At is approximately l/30th of the orbital time of the 
innermost planet. In order to verify the accuracy of the integra- 
tions, the maximum change in energy, e, permitted was 10" 4 . 
We define e as 

_ max\Ei -Eq| 
Eq 

where £, represents an individual measurement of the energy 
during the simulation, and Eq is the energy at time 0. There 
are two reasons for using this threshold in e. First, as the in- 
tegration scheme is symplectic, no long term secular changes 
will occur, so high precision is not required. Second, the sim- 
ulations needed to be completed in a timely manner. If a sim- 
ulation did not meet this energy conservation criterion, it was 
rerun with the timestep reduced by a factor of 10. The minimum 
timestep we used was P in „ er /3000. Despite this small timestep, 
a few simulations did not conserve energy, and were discarded, 
except that they were incorporated into the errors. Errors and 
error bars include information from unconserved simulations. 
Simulations which fail to conserve energy would most likely 
be labeled as unstable, as the failure of energy conservation un- 
doubtedly results from a close encounter between two bodies, 
which usually results in an ejection. Therefore the estimates 
for stability in the systems presented here should be considered 
upper limits. 

Throughout this paper we adopt the nomenclature of the dis- 
covery papers (planets have been labeled b, c, d, etc, with order 
in the alphabet corresponding to order of discovery). We will 
also introduce another scheme based on mass. Planets will be 
subscripted with a 1, 2, 3, etc, in order of descending mass. 
This new scheme is more useful in discussing the dynamics of 
the system. 

The short term simulations are integrated until one of the 
following criteria is met. 1) The simulation ejects a planet. 
Ejection is defined as the osculating eccentricity of one planet 
reaching, or exceeding unity. 2) The simulation integrates to 
completion at time r. For these simulations, r is defined as 

r = 2.8xl0 5 P!, (3) 

or 280,000 times the period of the most massive planet. This 
choice corresponds to 10 6 years for the v And system, as was 
simulated in Paper I. 

If a system ends without ejection, then the stability of the sys- 
tem must be determined. There are several possible definitions 
of stability. In Paper I, a system was stable if the osculating ec- 
centricity of each companion remained below 1 for the duration 
of the simulation. The most obvious flaw in this definition is 
that a planet could suffer a close approach and be thrown out to 
a bound orbit at some arbitrarily large distance. Such a system 
would bear no resemblance to the observed system, and hence 
should be labeled unstable. Here we adopt a more stringent def- 
inition, namely that the semi-major axes of all companions can- 
not change by more than a factor of 2. Changes in semi-major 
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axis represent a major perturbation to the system, therefore this 
second cut is conservative, and only eliminates systems which 
expel a planet to large distances without fully ejecting it. 

In addition to these short term simulations, we completed 
simulations to explore longer term stability (~ 10 8 yrs). For 
each system we ran ~ 10 simulations, chosen to cover a wide 
range of stable parameter space. For these runs we started 
with the final conditions of stable configurations and contin- 
ued them. These simulations therefore give us a handle on how 
the system is likely to evolve on timescales closer to its age 
(~ 10 9 yrs). Those systems which survived these longer runs 
are the best comparisons to the true system. Hence they are the 
best simulations for determining the factors that lead to plane- 
tary system stability. 

There are two notable problems with this methodology. First, 
we ignore the effects of general relativity, which may be impor- 
tant in some systems, specifically GJ876, and v And. General 
relativity was included in our treatment of v And in Paper I. In 
Paper I, the innermost planet had a negligible effect on the sys- 
tem, and we presume that general relativity will continue to be 
unimportant for the systems studied here. Second, we treat all 
particles as point masses. This is again especially troublesome 
for GJ876 and v And, because of their proximity to their (pre- 
sumably) oblate primaries. The sphericity of the stars also pro- 
hibits any tidal circularization of highly eccentric planets (Ra- 
sio et al. 1996). This may artificially maintain large planetary 
eccentricities, and increase the frequency of close encounters. 
However the eccentricities must become very large for this ef- 
fect to become appreciable. We therefore assume that this phe- 
nomenon will not adversely affect our results. Ignoring these 
two issues should impact the results minimally while speeding 
up our simulations considerably. 

2.3. Stability Maps 

When analyzing the simulations, it is useful to visualize the 
results in a stability map. In general a stability map is a 3 di- 
mensional representation of stability as a function of 2 parame- 
ters. In resonant systems, we find that several parameters deter- 
mine stability. The most important is the ratios of the periods 
of the two resonant planets, which we will call R: 

R=^-. (4) 

Pinner 

In coupled systems, e\ and e2 are the relevant parameters. The 
advantage to this visualization is that boundaries between sta- 
bility and instability are easily identified. 

The disadvantage of this form of visualization is that if the 
range of parameter space is not uniformly sampled (as it is 
here), we cannot visualize of the errors. It is therefore impor- 
tant to bear this disadvantage in mind. At the edges of stability 
maps, the data are poorly sampled and the information at the 
edges should largely be ignored. To aid the visualization we 
have smoothed the maps. If a bin contained no data then it 
was given the weighted mean of all adjacent bins, including di- 
agonal bins. This methodology can produce some misleading 
features in the stability maps. Most notable are tall spikes, or 
deep depressions in sparsely sampled regions. We comment on 
these types of errors where appropriate. 

The procedure as outlined overestimates the size of stable re- 
gions in two important ways. First, the integration times are 
generally less than 0.1% of the systems' true ages. As is shown 

2 http://obswww.unige.ch/~udry/planet/hd82943syst.html 



throughout this paper, instability can arise at any timescale. 
Therefore, the stability zone will continue to shrink as the sys- 
tem evolves. Second we have chosen a very generous cut in 
semi-major axis space. Other studies permit Aa to be no larger 
than 10% (Chiang, Tabachnik, & Tremaine 2001). Lowering 
this variation would undoubtedly also constrict stability zones. 

3. RESONANT SYSTEMS 

Two systems with orbital periods in 2:1 mean motion res- 
onance have been detected: HD82943 2 and GJ876 (Marcy et 
al. 2001a). Table 1 lists the orbital elements and errors for the 
resonant systems. For now we do not examine the 3:1 55Cnc 
system. The current best Keplerian fit to the observations put 
HD82943 and GJ876 just beyond perfect resonance. These 
planets all occupy high eccentricity orbits and hence have wide 
resonance zones. Simulations of these systems show that sta- 
bility is highly correlated with the ratio of the periods, R, and to 
a lesser degree on e\. These systems are the least stable as less 
than 20% of simulations survived to r. 

3.1. HD82943 

Two planets orbit the 1 .05 ± 0.05M© GO star (Santos et 
al. 2000) HD82943 at semi-major axis distances of 0.73 and 
1.16AU. Planet b is the inner and less massive, c, the outer and 
more massive. Gozdziewski & Maciejewski (2001) also exam- 
ined this system. They found that the system is most likely to 
be stable in perfect resonance. These simulations must run for 
340,830 years, t hd ^2943- The long term simulations were run 
for 100 million years. 

For this system we find that 18.8% ±4.3% of the trials were 
unperturbed to thd%29AZ, and 19.3% survived for 10 6 years, and 
4.5% failed to conserve energy. Fig. 1 shows the instability rate 
to 10 6 years. Most unstable simulations break our stability cri- 
terion within 10 4 years, but others survived more than 900,000 
years before being ejected. The asymptotic falloff to 10 6 years 
implies we have found most unstable configurations. In 91.3% 
of the trials, planet b, the inner and less massive planet, was 
ejected/perturbed. In order to check the simulations, Fig. 2 plots 
the rate of survival as a function of energy conservation. From 
this figure it appears that our limit of 10" 4 is reasonable, as there 
appears to be no trend in stability as a function of energy con- 
servation. The spike in survivability at 10~ 8 corresponds to reg- 
ular orbits that were stable for our initial choice of Af . The lack 
of a trend with energy conservation (specifically survival prob- 
ability increasing with decreasing e), implies our cutoff value 
of e is stringent enough. 

Stability in this system is correlated with R, e c , AM and AC3, 
where AM is the difference in mean anomaly, and AtES is the 
difference in initial longitude of periastron. Slightly beyond 
perfect resonance is the preferred state for this system. This 
system also requires the eccentricity of planet c to remain be- 
low 0.4. These features are shown in Figs. 3 and 4. In these 
greyscale images, black represents unsampled regions, darkest 
grey marks regions in which no configurations survived, grades 
of stability are denoted by lighter shades of grey, and white is 
fully stable. As with most stability maps in this paper, the outer 
2-3 gridpoints should be ignored. In Fig. 3, the R-e c stability 
map, the large "plateau" at low e, is therefore poorly sampled, 
as is the island at R=2.\5, e e =0.52. The most striking feature of 
this figure is that the best fit to the system, the asterisk, places 
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Table 1 - Initial Conditions for Resonant Systems 



System 


Planet 


Mass(Mj) 


Period(d) 


Eccentricity 


Long. ofPeri.(Deg) 


Time of Peri. (JD) 


GJ876 


c 


0.56 


30.12 ±0.02 


0.27 ±0.04 


330.0 ± 12.0 


2450031.4 ± 1.2 




b 


1.89 


61.02 ±0.03 


0.10 ±0.02 


333.0 ± 12.0 


2450045.2 ± 1.9 


HD82943 


b 


0.88 


221.6 ±2.7 


0.54 ±0.05 


138 ± 13 


2451630.9 ±5.9 




c 


1.63 


444.6 ± 8.8 


0.41 ±0.08 


96 ±7 


2451620.3 ± 12.0 
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it adjacent to stability. If R is changed by less that 5% the sys- 
tem has no chance of surviving even 1 million years. This map 
shows that the current fit to the system is not correct. However 
the elements do not need to change by much (specifically e c 
needs to be slightly lower) for the system to have a chance at 
stability. 
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Fig. 1. — The distribution of instability times for unstable configu- 
rations of HD82943. Most unstable systems survive for just 10 2 -10 4 
years before perturbations change a semi-major axis by more than a 
factor of 2. 
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Fig. 3.— The R-e c stability map for HD82943. The asterisk rep- 
resents the best fit to the system as of 31 July 2002. The data are 
most accurate closer to the asterisk. The system shows a clear bound- 
ary in phase space between unstable (dark grey) regions and stable 
(white) regions. Black represents unsampled data. The stable region 
at R = 2. 15, e c = 0.52 is a bin in which 1 out of 1 trials survived. 



The system also shows dependence on mean anomaly and 
longitude of periastron. Because of the symmetry of ellipses 
we will introduce a new variable, A, defined as: 
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Fig. 2. — Survival rate as a function of energy conservation for 
HD82943. The lack of a trend implies the results for the system are 
accurate. 



. _ J |aJi-aj 2 |, A<tt 

A -\ 360-|aJi-GJ 2 |, A>tt Pj 

where the subscripts merely represent 2 different planets, b and 
c for this system. The order is unimportant, as the we are only 
concerned with the magnitude of this angle. In Fig. 4, the A- 
AM stability map is presented. The asterisk marks the best fit 
to the system. Stability seems to follow a line represented by 

AM=1a+120=^(A + |). (6) 

This relation is purely empirical. As is shown in the following 
sections, this interdependency is unusual for extra-solar plan- 
etary systems. Similar plots of AM or A versus R, show the 
same R dependence as in Fig. 3. Therefore R is clearly the most 
important parameter in this system, but these other 3 also play 
an important role in the system. As more observations of the 
system are made, HD82943 should fall into the region defined 
by 1 .95 < R < 2.1, e c < 0.4, and Eq. 6. 
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becomes locked at A w 100°. Although not shown, after these 
initial perturbations the system settles down to a configuration 
in which a b « Q344AU, e b « 0.79 and a c w 2400AU, e c w 0.99. 
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Fig. 4.— The A-AM stability map for HD82943. The asterisk rep- 
resents the best fit to the system as of 31 July 2002. The data are most 
accurate closer to the asterisk. Stability appears to follow the line rep- 
resented by Eq 6. Note however that the system is also constrained to 
A < 75° and AM > 30°. The island at A = 80°, AM = 20° is a bin in 
which 1 out of 1 trials survived. 

HD82943 shows a wide range of dynamics. Some examples 
of these are shown in Figs. 5-8. The initial conditions of these 
systems are presented in Table 2. Fig. 5 shows an example of 
the evolution of a regular system, HD82943-348, which shows 
no evidence of chaos. Note that instead of A(t), we present the 
distribution function of A, P(A), the probability of A, vs. A. 
This representation of A shows that the motion is like that of 
a harmonic oscillator; the longitudes of periastron are librating 
with an amplitude of 60°. 

Fig. 6 (HD82943-382) is a stable case which is clearly 
chaotic. Although the eccentricities remain close to their initial 
values, the inclinations jump to large values quickly. Note that 
A never exceeds 75°, but its motion is slightly nonharmonic, 
another indication of chaos. 

In Fig. 7 we plot the evolution of a system which ejects planet 
b. This system experiences close approaches very quickly as is 
evidenced in the top left panel of Fig. 7. The eccentricities un- 
dergo dramatic fluctuations, and the inclinations nearly reach 
40°. Although poorly sampled, A suggests circulation. 

In Fig. 8, the orbital evolution of simulation HD82943-216 
is shown. This system perturbed planet c after 280,000 years, 
and despite the high eccentricities the system obtained (>0.75 
for both planets), remained bound for 10 6 years. The inclina- 
tions also show large growth. Although initially A = 41°, the 
system becomes locked at just over 100°. This is misleading 
as this is the A distribution for the full 10 6 years. It is actu- 
ally the superposition of 2 modes. Initially, until perturbation at 
280,000 years, the system circulates with a slight preference to- 
ward anti-alignment. However after the perturbation the system 
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Fig. 5.— Orbital evolution of HD82943-348, a stable, regular con- 
figuration. The data here are smoothed over a 20,000 year interval. 
Top Left: The evolution in semi-major axis. The planets show slight 
variations due to resonant interactions. Top Right: The eccentric- 
ities oscillate with a 700 year period, with 0.62 < e b < 0.85 and 
0.05 < e c < 0.45. This short timescale is not visible in this plot. 
Bottom Left: The inclinations experience oscillations on 2100 year 
timescale, again too short to be visible in this plot. The ranges of 
inclination are 1° < i;, < 7° and < i c < 3°. Bottom Right: From 
this distribution function we see the lines of node librate harmoni- 
cally with an amplitude of 60°. The dashed vertical line represents 
Ao, the initial value of A. 
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Table 2 - Selected Simulations of HD82943 



Trial 


eb,o 


e c ,o 


Ro 


Ao (°) 


e 


Comments" 


216 


0.563 


0.345 


2.0108 


42.2 


2.6 x 10" 7 


C,P(c,299.9) 


348 


0.617 


0.419 


1.9656 


13.6 


1.5 x 10~ 8 


R 


382 


0.516 


0.322 


2.0303 


32.4 


9.2 x 10~ u 


C 


698 


0.481 


0.440 


2.0948 


68.4 


1.3 x lO" 10 


C,P(b,7.5),E(b,10.5) 



a R=Regular, C=Chaotic, P=Perturbed(Planet, Time(10 3 yr)), E=Ejected(Planet, 
Time(10 3 yr)) 
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Fig. 6.— Orbital evolution of HD82943-382, a stable, chaotic con- 
figuration of HD82943. The data are averaged on a 10,000 year in- 
terval. Top Left: The evolution in Semi-Major Axis. The planets 
clearly never experience a close encounter. Top Right: The eccentric- 
ities experience low amplitude (w 13%) chaotic oscillations. Bottom 
Left: The inclinations initially jump to large values and experience 
large amplitude (« 70%) fluctuations. Bottom Right: The A distribu- 
tion function suggests that A is generally librating, but that there are 
chaotic fluctuations superimposed on this motion. The dashed vertical 
line is the value of Ao. 
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Fig. 7.— Orbital evolution of HD82943-698, the perturbation and 
ejection of HD82943b. The orbital parameters are sampled every 100 
years. Top Left: The planets experience close encounters immediately 
as changes in semi-major axis are visible within 500 years. Top Right: 
The eccentricities experience large amplitude fluctuations, culminat- 
ing in the ejection of planet b just after 10,000 years. Bottom Left: 
The inclinations vary slightly until a close encounter at 1100 years 
which sends both inclinations up substantially. The planets remain in 
this state until just prior to ejection, when they return to coplanarity. 
Bottom Right: The A distribution function suggests that A remains 
below 100°, but may occasionally circulate. 
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Fig. 8.— Orbital evolution of HD82943-216, the perturbation of 
HD82943c. The parameters are averaged on 10,000 year intervals. 
Top Left: The planets' semi-major axes are stable and show no signs 
of close encounters until 260,000 years. At this point planet c actually 
crosses b's orbit. This initial encounter leads to more encounters as 
a c reaches 3AU by 280,000 years, tripping the criterion for instability. 
Top Right: The eccentricities experience secular change until 210,000 
years. The system then moves into a lower eccentricity state. The 
eccentricity then grows to large values and remain at their final val- 
ues for another 700,000 years. Bottom Left: As with eccentricity the 
inclinations show slow secular change until 210,000 years. The incli- 
nations then leap up to 30° in the case of planet b. Bottom Right: The 
A distribution function is the sum of 2 motions: the pre-perturbation 
motion is circulation, the post-perturbation motion is fixed close to 
110°. The dashed line represent Ao. 
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We ran 10 simulations for 10 years. The initial conditions 
and results of these simulations are presented in Table 3. The 
eccentricity evolution of 4 simulations is shown in Fig. 9. These 
simulations, which ran for 100 million years, show some con- 
figurations are regular (top left), some are chaotic (top right, 
bottom left), and that instability can arise at any timescale (bot- 
tom right). These long term simulations show that regions exist 
in phase space in which this system can survive for billions of 
years. 



"o 20 40 60 80 20 40 60 80 100 

Time (10 s Yr) Time (10 6 Yr) 

Fig. 9. — Long term simulations of the HD82943 system. These 
data are averages over 50,000 year intervals. Top Left: Evolution of 
HD82943-035. An example in which the system is regular. Top Right: 
Evolution of HD82943-021, an example of chaotic evolution. Bottom 
Right: Evolution of HD82943-000, another example of chaotic mo- 
tion. Bottom Left: Evolution of HD82943-032. A chaotic system 
which ejects planet b after 2.4 million years (7t#d82943). 
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Table 3 - Results of Long Term (10 8 yr) Simulations of HD82943 



Trial 


eb,o 


e c ,o 


Ro 


Ao (°) 


e 


Result" 


000 


0.561 


0.395 


2.0269 


26.5 


5.2 x 10~ 8 


C 


021 


0.571 


0.333 


2.0615 


39.0 


1.0 x 10~ 8 


C 


032 


0.414 


0.556 


2.0556 


33.3 


4.95 


e(2.4),E(b,2.4) 


035 


0.511 


0.236 


2.0322 


6.7 


6.8 x 10- 9 


R 


036 


0.534 


0.197 


2.0080 


29.5 


6.1 x 10-" 


R 


040 


0.600 


0.393 


2.0766 


52.8 


0.22 


e(1.005),E(c,1.14) 


043 


0.542 


0.386 


2.0223 


24.8 


2.4 x 10~ 8 


C 


057 


0.549 


0.348 


2.0382 


27.0 


1.1 x 10" 8 


C 


073 


0.468 


0.265 


2.0626 


56.8 


6.3 x 10- 8 


R 



a R=Regular, C=Chaotic, e=Energy conservation failed (Time (10 6 yr)), E=Ejection(Planet, 
Time (10 6 yr)) 
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The Stability of Planetary Systems 



Several papers have suggested that secular resonance locking 
maintains stability in ESPS with large eccentricities (Rivera & 
Lissauer 2000, Rivera & Lissauer 2001a, Chiang, Tabachnik, & 
Tremaine 2001). Specifically they suggest that the orientation 
of the planets' ellipses should be aligned (A «0). By examin- 
ing A in stable, regular long term simulations we can determine 
if the longitudes of periastron remain locked. The probability 
distribution of A for these same four long term simulations is 
presented in Fig. 10. A is sampled once every 100 years. The 
top left plot of Fig. 10 is similar to that of an harmonic oscilla- 
tor; this configuration is librating about A=0 with an amplitude 
of 40°. The other plots are systems which are not librating, but 
instead show more random motion. From this figure we see that 
regular motion is correlated with libration about alignment, but 
chaotic and unstable generally shows chaotic A behavior. 



of parameter space is stable. On 10 6 year timescales, 2.4% of 
configurations survived, but only 1 .7% were still unperturbed, 
and 5% failed to conserve energy. In unstable cases planet c, 
the inner and less massive, was ejected/perturbed nearly 96% 
of the time. Fig. 1 1 shows the instability rate. GJ876 shows the 
same trend as HD82943; most unstable configurations break up 
in just a few hundred dynamical times. Again the asymptotic 
nature of this plot implies most unstable configurations have 
been identified. This system shows the same sort of trend in e 
as HD82943, which proves that our results our correct. 



50 100 150 
A(°) 



50 100 150 
A(°) 



Fig. 10. — The distribution function of A (sampled every 100 years) 
for four stable cases of HD82943. These four systems are the same as 
in Fig. 9. The top left plots a system librating about A=0. The other 
plots show that chaotic motion is usually associated with a circulating 

A. 




log 10 (t 



Fig. 11. — The ejection rate for GJ876. In this system unstable 
configurations are usually ejected within 100 years. The rate asymp- 
totically falls to zero by 10 5 years. 



3.2. GJ876 

Two planets orbit the 0.32±0.05M Q (Marcy et al. 2001a) 
M4 star GJ876, also known as Gliese 876. This system is very 
similar to HD82943, the major difference being that the plan- 
ets lie closer to their primary. The semi-major axes of these two 
planets are 0. 13 and 0.21 AU. Note that in this system planet c is 
the inner and less massive planet. A substantial amount of work 
has already been done on this system. Notably, in the discovery 
paper (Marcy et al. 2001) that stable configurations exist in the 
system. Rivera & Lissauer (2001b) show that Keplerian fitting 
for this system is not precise enough to accurately determine 
the orbits. They suggest, through N-body fitting, that GJ876 
must actually lie in perfect resonance, and that the orbital el- 
ements provided in the discovery paper (which are used here) 
suffer from a systematic error. Because Pi is so short for this 
system (60 days) the evolution of the orbital elements has been 
observed, and corroborated this theory. This section therefore 
serves as a check to this hypothesis. 

We ran GJ876 for 10 6 years, but t C j%k, corresponds to only 
47,000 years. On 47,000 year timescales, only 5.6% ± 2.8% 



Unlike HD82943, there are no obvious zones of stability. We 
only see isolated islands in the R-et stability map presented 
in Fig. 12. We choose these parameters for our stability map as 
they were the strongest indicators of stability in HD82943, and 
because of the suggestion that the system actually lie in per- 
fect resonance. Close to R = 2.00 we sampled two simulations 
near R = 2.02 and e/,=0.7. Both these were stable. However 
with such poor statistics, and at such a large (relative) distance 
from perfect resonance, we cannot comment on the likelihood 
that the system would be more stable in perfect resonance. We 
can, however, point out that there are isolated configurations 
which may hold stable orbits, and prolonged observations of 
this system should demonstrate whether it is indeed in perfect 
resonance. However this lack of a large stable region strength- 
ens the hypothesis that this system lies in perfect resonance. 
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Fig. 13. — The dependence of stability on initial A. The data above 
50° are poorly sampled, but with such a large difference between their 
values and the mean of 5.7%, they do suggest that stability may be im- 
proved at A > 50". 
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Fig. 12. — TheR-eb stability map for GJ876. The asterisk marks the 
best fit to the system as of 7 Aug 2002, and the values for stability are 
more accurate closer to the asterisk. In this system, as in HD82943, 
the two relevant orbital elements are ei and R. There are no contigu- 
ous regions of stability, only small isolated pockets which may hold 
stable zones. 

This system does, however, lie very close to A=0. However 
this proximity to alignment has no bearing on the stability of the 
system, in fact, it may actually diminish its chances of stability. 
In Fig. 13 the probability of stability as a function of initial A 
is shown. Although the values for large A are poorly sampled, 
four data points at 100% stability does suggest larger values of 
A may be more stable. 

Similar dynamics are present in GJ876 as in HD82943. In 
Figs. 14-17, we show 4 examples of stable and unstable config- 
urations. The initial conditions for these sample simulations are 
listed in Table 4. In Fig. 14 the orbital evolution of simulation 
GJ876-904 is shown. This system is one of very few (<10) sim- 
ulations which show approximately regular motion. The semi- 
major axes show resonant perturbations, and A shows a motion 
typical of chaos; it appears to be the superposition of libration 
and circulation. 

In Fig. 15 the evolution of a stable, yet chaotic, configura- 
tion is plotted. Although this system was stable for 47,000 
years, planet b was perturbed after 220,000 years. This system 
however remained bound for 10 6 years. The large eccentric- 
ity oscillations continue on to 10 6 years, and planet c tends to 
remain in a retrograde orbit. After 10 6 years a c = 0.05 15AU, 
0.06 < e c < 0.85, i c « 120°, a b w 0.87AU, e b « 0.77, and 
h < 5°. The period of e c oscillations remains at 3300 years. 
The A evolution further belies the chaos in this system, as it 
tends toward anti-alignment, but also circulates. As before this 
distribution is over 10 6 years, but as the system remains in ap- 
proximately the same state from 25,000 years to 10 6 years, this 
plot is a fair representation of the motion during the first Tame- 
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Fig. 14. — Orbital evolution of GJ876-904, the stable, regular mo- 
tion of GJ876. The data are smoothed on 3000 year intervals, which 
eliminates all short period behavior. Top Left: The low order res- 
onance and large eccentricities induce low amplitude oscillations in 
the semi-major axes of the planets with a period of 2 years. This short 
frequency is not visible in this plot due to the smoothing interval. Top 
Right: The evolution of eccentricity for this configuration is that of 
two eigenmodes, also on a period of 2 years. Bottom Left: The in- 
clinations vary slightly with a frequency of 250 years. Bottom Right: 
The A distribution function suggests that A sometimes librates with 
an amplitude of » 80°, and sometimes circulates. The dashed line 
marks Ao. 



The Stability of Planetary Systems 



Table 4 - Selected Simulations of GJ876 



Trial 


e c ,o 


eb,o 


Ro 


Ao (°) 


e 


Comments 


029 


0.177 


0.073 


2.0241 


11.6 


4.6 x 10" 5 


C 


563 


0.239 


0.122 


2.0247 


3.65 


2.2 x 10~ 8 


C, P(b,5.1), E(b,5.3) 


858 


0.338 


0.107 


2.0255 


13.1 


6.2 x 10" e 


C, P(c,14.9), E(c,151.9) 


904 


0.351 


0.072 


2.0303 


51.0 


2.7 x 10~ 9 


R? 



° R=Regular, C=Chaotic, P=Perturbed(Planet, Time (10 3 yr)), E=Ejected(Planet, 
Time (10 3 yr)) 



1 
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Fig. 15. — Orbital evolution of GJ876-029, a chaotic stable config- 
uration of GJ876. Top Left: Although the semi-major axes vary, they 
don't change by a factor of 2 until 220,000 years=4.7roj876. Top Right: 
The remarkable eccentricity evolution of this system. These oscilla- 
tions persist for 10 6 years. Bottom Left: The inclination of planet 
b varies a small amount, generally staying below 5°. Planet c how- 
ever experiences wild fluctuations, however it does eventually settle 
to 120". Bottom Right: This curious A distribution function suggests 
that A prefers anti-alignment. This implies that a protection mecha- 
nism is keeping the system from breaking apart despite the extremely 
large values of et- 



In Fig. 16, the evolution of GJ876-563 is shown. This con- 
figuration ejects planet c in just over 5000 years. The system 
shows evidence for close approaches right from the beginning 
as the semi-major axes show small amplitude, chaotic fluctua- 
tions. The inclinations and eccentricities also experience a large 
degree of chaos. Although initially A is very close to 0, it cir- 
culates for the duration of the simulation. 




Fig. 16.— Orbital evolution of GJ876-563, the ejection of GJ876c. 
Top Left: The semi-major axes undergo small fluctuations until a vio- 
lent encounter ejects the inner planet after 5000 years. Top Right: The 
eccentricities experience large amplitude fluctuations, until e c eventu- 
ally reaches unity. Bottom Left: The inclination of planet b initially 
jumps and remains at over 30°, but then returns to coplanarity just 
before ejection. This type of motion was also seen in Fig. 8. Bottom 
Right: As with the majority of chaotic configurations, A circulates. 
Although poorly sampled, A appears to have no preferred value. 



In Fig. 17 the evolution of a system which perturbs the outer 
planet in just 9000 years is shown. The system ejects planet b in 
152,000 years. Before reaching tg/876 this system experiences 
some remarkable evolution in a, e, and i. Note, too, that A very 
quickly moved to an anti-aligned configuration. 
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Fig. 17.— Orbital evolution of GJ876-858, the perturbation of 
GJ876b. Planet b was eventually ejected after 152,000 years. Top 
Left: The semi-major axes evolve quiescently for 9000 years, until a 
close approach increases at, marked by the p. Although at returns 
to its initial value by tgjkk, just prior to ejection at reached 750AU. 
Top Right: The eccentricities immediately jump to very large values. 
e c varies wildly between 0.6 and 0.99. After b is kicked out to large 
a, the oscillations become much smaller. For nearly 10,000 years et 
remains above 0.98, but it doesn't reach unity until 152,000 years. 
Bottom Left: As with eccentricity, the inclination of c jumps wildly 
for 9000 years, even reaching 162° just prior to perturbation. How- 
ever it is unclear if this large inclination produces the perturbation. 
Bottom Right: As with the e and (', A immediately moves from its 
starting position. However A remains very close to an ti- alignment for 
the duration of the simulation. 



Long term simulations for this system were integrated for 
27.5 million years. A complete summary of the long term sim- 
ulations for this system is presented in Table 5. Fig. 18 plots the 
eccentricity evolution of 4 simulations. Some systems appear 
regular throughout (top left). Some configurations are chaotic 
for the duration of the simulation (top right, bottom left), and 
others may eject a planet after an arbitrarily long period of time 
(bottom right). 
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Table 5 - Results of Long Term (27.5 x 10 6 yr) Simulations of GJ876 



Trial 


e c ,o 


eb,o 


Ro 


Ao (°) 


e 


Result" 


290 


0.0884 


0.248 


2.02353 


45.3 


2.4 x 10~ 9 


C 


300 


0.108 


0.161 


2.02641 


15.2 


0.5 


e(4.06), E(b,4.08) 


362 


0.124 


0.366 


2.02624 


29.3 


2.6 x 10~ 8 


R 


404 


0.121 


0.246 


2.02548 


46.43 


0.6 


e(17.6), C 


609 


0.0875 


0.241 


2.02342 


53.5 


2.5 x 10~ u 


R 


655 


0.110 


0.256 


2.02601 


13.9 


0.19 


e(11.2), P(b,24.0) 


661 


0.105 


0.208 


2.02778 


33.1 


0.18 


e(2.4), E(c,2.4) 


739 


0.107 


0.273 


2.02726 


31.2 


0.02 


e(1.4), E(b,1.4) 


895 


0.0816 


0.228 


2.02242 


19.4 


2.3 x 10- 8 


C 


948 


0.0628 


0.204 


2.02803 


4.80 


4.1 x 10" 4 


e(0.25), C 



a R=Regular, C=Chaotic, e=Energy conservation failed (Time 10 6 yr)), P=Perturbed 
(Planet, Time (10 6 yr)), E=Ejected (Planet, (10 6 yr)) 
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Time (10 8 Yr) Time (10 fi Yr) 

Fig. 18. — Eccentricity of 4 long term simulations of GJ876. Top 
Left: Evolution of GJ876-362, a regular configuration. Top Right: 
Evolution of GJ876-290, a chaotic, yet stable configuration. Bottom 
Left: Evolution of GJ876-895, a chaotic, stable configuration. Bottom 
Right: Evolution of GJ876-300 which ejects planet b after 4 million 
years. 

Fig. 19 shows the distribution function of A for the same four 
systems. The regular system (top left) shows a configuration 
which usually librates with an amplitude of 80°, but with oc- 
casional circulation. The two chaotic examples (top right, bot- 
tom left) have flat distribution functions. Not surprisingly the 
unstable trial shows a very peculiar distribution function. As 
in HD82943 we see that regular systems tend to librate, and 
chaotic configurations circulate. 
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Fig. 19. — The A evolution of the same four simulations in Fig. 19. 
Chaotic systems (top right, bottom left) show no signs of libration, 
while regular systems (top left) are librating, but with occasional cir- 
culation. The unstable example (bottom right) shows a very strange 
distribution, which only demonstrates the chaos of this system. 
3 http://ssd.jpl. nasa.gov/elem_planets. html 



Although the evidence is compelling that GJ876 does in fact 
lie in perfect resonance, our work demonstrates that stable, reg- 
ular systems do exist close to the observed Keplerian fit. More 
observations of this system will demonstrate if the system is in 
perfect resonance. This work clearly demonstrates that stable 
regions do exist for a system like GJ876 just beyond perfect 
resonance. Should this system lie in perfect resonance, then 
this work shows that unstable regions lie very close to its con- 
figuration. 

New astrometric data has confirmed the mass of the outer 
planet in this system (Benedict et al. 2002). This therefore pro- 
vides the only system with a known mass. The plane of b's orbit 
is inclined by 6° to the line of sight. Benedict et al. confirm the 
mass and semi-major axis of this planet to be statistically iden- 
tical to those presented in Table 1 . However for this paper, the 
lack of data for planet c precludes any new insights into the dy- 
namics of the system. At best, if the system is approximately 
coplanar, then the variations used here are indeed representative 
of the true system, and our results are more robust. 

4. INTERACTING SYSTEMS 

Four known systems meet the interacting system criterion: v 
And (Butler et al. 2000), 47UMa (Fischer et al. 2002), the SS 3 
and HD12661 (Fischer et al. 2003). In this paper we will limit 
ourselves to the first three. The number, placement, and size 
of the planets in each system are quite different, but all have 
at least two planets that lie in between the 2:1 and 10:1 reso- 
nances, v And was the first known ESPS, and was the subject 
of Paper I. The experiment in Paper I is the procedure for this 
paper, and the simulations have been performed again. 47UMa 
was announced in 2001, and, at first glance, appears more like 
the SS than v And. Performing this experiment on the SS is 
problematic. The errors in the orbital elements of the SS are 
drastically smaller than for the ESPS, therefore fitting the SS 
into the procedure requires inflating the SS orbital element er- 
rors to values comparable to those of the ESPS. Essentially we 
are asking what would we observe if we took radial velocity 
measurements of our sun. We will compare the ESPS to both 
the gas giant system (§4.3) and the Jupiter-Saturn system (§4.4). 
These three coupled systems' orbital elements are summarized 
in Table 6. Interacting systems show broad regions of stability 
which are correlated with eccentricity. 

4.1. v Andromedae 

The v And system is a combination of a separated system and 
an interacting system. Three planets orbit the 1.02±0.03M Q 
(Gonzalez & Laws 2000) F8 star v And. The inner planet, b 
orbits at 0.04AU. The other planets, c and d, orbit at larger dis- 
tances (0.8AU and 2.5AU respectively), but are significantly 
more eccentric. The outer planet is the most massive, therefore 
TvAnd corresponds to 10 6 years. 

The v And system was the subject of Paper I, and has been 
the focus of intense research since its discovery. The apparent 
alignment of the apses of planets c and d has sparked the most 
interest with several groups claiming that the system must be 
secularly locked, or at least librate about A=0 (Rivera & Lis- 
sauer 2000; Rivera & Lissauer 2001b; Chiang, Tabachnik, & 
Tremaine 2002), while others suggest this alignment may be 
a chance occurrence (Paper I; Stepinski, Malhotra, & Black 
2000). However, these groups and others (Laughlin & Adams 
1999) all agree that this system, as presented, can be stable for 
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Table 6 - Initial Conditions for Interacting Systems 



System 


Planet 


Mass(Mj) 


Period(d) 


Eccentricity 


Long, of Peri.(Deg) 


Time of Peri. (JD) 


SS 


Jupiter 


1.000 


4331.6 ±43.3 


0.0484 ±0.1 


14.8 ±90.0 


2449896.3 ±25.0 




Saturn 


0.297 


10759.7 ± 107.8 


0.0542 ±0.1 


92.4 ±90.0 


2450411.1 ±25.0 




Uranus 


0.0459 


30704.9 ±307.0 


0.0472 ±0.1 


171.0 ±90.0 


2447230.0 ±25.0 




Neptune 


0.0541 


60197.2 ±602.0 


0.0086 ±0.1 


45.0 ±90.0 


2442071.3 ±25.0 


v And 


b 


0.69 


4.61706 ±0.0003 


0.015 ±0.015 


32.0 ± 243.0 


2459991.588 ±3.1 




c 


1.96 


241.1 ±1.1 


0.25 ±0.11 


251.0 ±33.0 


2450160.1 ±20.8 




d 


3.98 


1309.0 ±30.0 


0.34 ±0.11 


255.0 ± 17.0 


2450044.0 ±40.5 


47UMa 


b 


2.54 


1089.0 ±3.0 


0.06 ±0.014 


172.0 ± 15.0 


2453622.0 ±34.0 




c 


0.76 


2594.0 ±90.0 


0.1 ±0.1 


127.0 ±56.0 


2451363.5 ±493.0 
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at least 10 years. 

On a 10 6 year timescale, 86.1%±3.3% of simulations sur- 
vived, and 0.4% failed to conserve energy. This value is less 
than la from the value published in Paper I, 84.0%±3.4%. 
Fig. 20 shows the perturbation rate as a function of time. Once 
again we see that most unstable configurations eject a planet im- 
mediately, and the rate falls to 4% by 10 6 years. The fact that 
ejections occur right up to 10 6 years implies that we have not 
detected all unstable situations, and that the stability map for 
this system contains more unstable configurations, and hence 
the plateau is smaller and/or the edge is steeper after 10 9 years. 

There is one notable difference between the results of Paper 
I, and those reported here: the frequency of ejections of each 
planet is different. In Paper I planet b was ejected 10% of the 
time, c, 60%, and d, 30%. The SWIFT runs ejected planet b 
39% of the time, c, 54%, and d, 7%. Fig. 21 plots the survival 
probability in this system as a function of energy conservation. 
Stability peaks at e = 10~ 8 but quickly drops. Although this plot 
is qualitatively different from Fig. 2, we again note that this im- 
plies the simulations are valid. This plot is typical for the in- 
teracting systems, confirming our hypothesis that we need only 
maintain e < 10" 4 for the duration of every simulation. 
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Fig. 20. — The ejection rate of unstable configurations in v And. 
Configurations may eject planets right up to 10 6 years. It is therefore 
unclear how many more configurations might become unstable. 

In Paper I, a stability map in e c and was presented in Table 
2. This table showed that and, to a lesser degree e c deter- 
mined the stability of the system. In Fig. 22, the v And stability 
map is presented. It is nearly an exact match to that in Paper 
I. However the best fit to the system 4 has changed since then 
and the system has moved away from the edge slightly. It is 
important to note that two different integration schemes pro- 
duced the same results. From Fig. 22, it is clear that v And 
lies near instability. The edge of stability in v And, however, is 
fundamentally different than in resonant systems. In this inter- 
acting system, a large region of phase space is fully stable, the 
"plateau", but there is a sharp boundary in eccentricity space, 
the "edge", beyond which the system quickly moves into a fully 
unstable region, the "abyss". Although it appears that both mor- 
phologies are on the edge of stability, they are different types of 

4 http://exoplanets.org/upsandb.html 
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Fig. 21. — Survival probability as a function of energy conservation 
in v And. The instability at low e implies the results for v And are 
robust. 
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Fig. 22. — Stability map for v And. Stability in this system depends 
on ej and, to a smaller degree, on e c . As in resonant system stability 
maps, precision is correlated with distance from the asterisk, which 
marks the best fit to the system as of 24 Sept 2002. v And lies near 
the edge of stability. 

As previously mentioned there are some unusual features of 
this system; the lines of node are nearly aligned, and the system 
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lies close to the 11:2 mean motion resonance. This is a weak 
perturbation, but between these 2 massive planets, this may be 
important. However a quick inspection of plots of stability ver- 
sus A, Fig. 23, and R, Fig. 24, shows that there is no statistically 
significant affect caused by these to (potential) resonances. We 
do note that our integration time may not be long enough to 
detect the importance of the 11:2 resonance. 



actually demonstrates the breakdown of our model. Planet b's 
eccentricity oscillates with an amplitude of 0.3 and a period of 
120,000 years. Unfortunately, v And b is tidally locked by its 
parent star with a period of 0.01 1 years. The timescale for tidal 
circularization is 8 x 10 7 yrs (Trilling 2000). We address this 
potential inconsistency in §7. A for this system librates with an 
amplitude of 50°, an indicator of regular motion. 




60 

A (o) 

Fig. 23. — Stability as a function of A in v And. Although the best fit 
to the system places it very close to alignment, marked by the dashed 
vertical line, there is no significant trend with A. 




Fig. 24. — Stability as a function of R in v And. There is no evidence 
that this weak mean motion resonance affects the stability of v And. 
Although there are fully stable points outside of resonance, all vary 
by less than la of the mean stability rate of 0.86. 

As has been shown in other papers, this system exhibits both 
regular, and chaotic motion. In Figs. 25-29, we present 5 ex- 
amples of this system. The initial conditions of these config- 
urations are listed in Table 7. In Fig. 25, the orbital evolution 
of a regular, stable configuration is plotted. However, this plot 




Fig. 25. — The orbital evolution of v And-076, a stable, regular con- 
figuration, smoothed on 25,000 year intervals. Top Left: The semi- 
major axes show no evidence of perturbations. Top Right: The eccen- 
tricities experience simple sinusoidal variations. The period of planet 
b oscillations is 100,000 years, while c and d oscillate on a 7000 year 
period. The large amplitude of e b is most likely unphysical due to 
tidal locking with the parent star. The apparently irregular behavior 
of e\, is an artifact of its long cycle. Bottom Left: The inclinations are 
also regular, although planet b's inclination is affected by both planets 
on its 20,000 year period. Planets c and d oscillate in inclination on a 
4000 year period. As in eccentricity, the slightly chaotic appearance 
of ib is an artifact of the sampling time convolved with the physical 
period. Bottom Right: The A distribution librates with an amplitude 
of 50°. 



In Fig. 26 the evolution of a stable, but chaotic system, v 
And-236, is presented. In this configuration ad is mildly per- 
turbed, but the eccentricities and inclinations undergo wild fluc- 
tuations. A also shows evidence of chaos as it librates and cir- 
culates. 
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Table 7 - Selected Simulations of v And 



Trial 


e c ,o 


ed,o 


Ro 


Ao (°) 


Comments" 


000 


0.276 


0.486 


5.583 


38.3 


C, P(c,218), E(b,218) 


014 


0.177 


0.470 


5.509 


20.2 


C, P(d,219), E(d,219) 


020 


0.299 


0.532 


5.207 


21.7 


C, P(d,10.1), E(c,31.3) 


076 


0.394 


0.375 


5.566 


11.1 


R 


236 


0.168 


0.428 


5.386 


81.0 


C 



° R=Regular, C=Chaotic, P=Perturbed(Planet, Time(10 3 yr)), E=Ejected(Planet, 
Time(10 3 yr)) 
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Fig. 27. — The orbital evolution of v And-000, the ejection of v 
And b. Top Left: The system appears to be experiencing perturbations 
right from the start as ad evolves chaotically. Eventually a close ap- 
proach between b and c produces an ejection. Top Right: Once again 
we see that tidal circularization would alter the outcome of this trial. 
Although it is unclear if circularization would have been enough to 
prevent the close approach that ejects planet b. Bottom Left: In con- 
trast to eccentricity, inclination evolves regularly for 200,000 years, 
when a close approach sends b into a retrograde orbit, while i c and id 
continue to evolve sinusoidally. Bottom Right: This unusual A dis- 
tribution function reveals further chaos in the system, and, as is seen 
in other systems, appears to be the sum of a libration mode, and a 
circulation mode. 

In Fig. 28 we show the orbital evolution of a system which 
perturbs planet d, but ejects planet c. The behavior of the eject- 
ing planet reaching very large semi-major axis distances, and 
subsequently returning, only to be ejected was also seen in Fig 
17. Note, too, that the peak in eb corresponds with 



Fig. 26. — The orbital evolution of v And-236, a stable, chaotic con- 
figuration. Top Left: The semi-major axes show perturbations of plan- 
ets c and d. Top Right: The eccentricities are quite chaotic, especially 
e b . However, we see very large eccentricities, and low semi-major 
axes for planet b. This again suggests that these large values may be 
unphysical. However, this does not change the fact that this system 
is chaotic. Bottom Left: The inclinations also show a large degree of 
chaos. As with the eccentricities, it is b which shows the most devia- 
tions from regular motion. Bottom Right: The A distribution suggests 
that this system librates and circulates. This is the typical behavior of 
a chaotic system. 



v And may eject any planet. However, as demonstrated 
in Fig. 25, the model breaks down for large e/,. The path to 
ejection for one such ejection of planet b is shown in Fig. 27. 
The eccentricities stay generally low (<0.2) for nearly 200,000 
years, but prior to ejection e b grows to 0.99 at times, before 
finally being ejected. 
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Fig. 28. — The orbital evolution of v And-020, the ejection of v And 
c. Top Left: The semi-major axes evolve quiescently for 10,000 years 
before perturbing planet d, marked by the p. 10,000 years later planet 
c is perturbed to over 120AU. The planet then returns to low a, but is 
quickly ejected after another encounter with planet d. Top Right: This 
configuration experiences wild oscillations from the very beginning. 
Bottom Left: The inclinations also suffer large, chaotic fluctuations. 
Shortly after planet d is perturbed, planet b experiences a short period 
of retrograde motion. Bottom Right: Although poorly sampled, this 
graph clearly shows that A evolves chaotically. 

Finally we show the ejection of planet d, the most massive 
planet in the system in Fig 29. This system evolves quasi- 
regularly for 200,000 years before suddenly ejecting planet d. 
The only hints of chaos are in the evolution of e c and e^, al- 
though there are some slight asymmetries in the sinusoidal evo- 
lution of e b . A is generally librating, although in the final 
20,000 years it does circulate as the system destabilizes. 

Long term simulations run for 100 million years. Fig. 30 is 
the eccentricity evolution of 4 simulations. In one case (bot- 
tom right panel of Fig. 30) the inner planet is ejected after 55 
million years. The top left panel shows a system undergoing 
chaotic evolution. The other two panels show regular motion. 
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Fig. 30. — The long-term eccentricity evolution of 4 simulations of 
v And. Top Left: Eccentricity evolution of v And-006. This is an 
example of regular evolution. Top Right: Eccentricity evolution of v 
And-054. This is a chaotic configuration. The chaos has been mostly 
smoothed over, though, as the data represent 10,000 year averages. 
Bottom Left: Eccentricity evolution of simulation v And-288. An- 
other chaotic configuration. Bottom Right: Eccentricity evolution of 
v And-192. A chaotic system which ejects the inner planet after 55 
million years. This evolution is suspect as the effects of tidal circular- 
ization undoubtedly play a role in the evolution of et- 

Fig. 31 shows the A distribution of these configurations. Sec- 
ular resonance locking does not occur in the v And system; the 
ellipses tend to be anti-aligned. This is the same result as in 
Paper I. However there do appear to be several different modes 
of stability for the apses in this system. The panels in Fig. 31 
correspond to those in Fig. 30, therefore the top left is the regu- 
lar case. The chaotic systems (top right, bottom left) show the 
A distribution signature of chaos, as does the example which 
ejects planet b (bottom right). In Table 8 we present a summary 
of all long term simulations for v And. 



Fig. 29. — The orbital evolution of v And-014, the ejection of v 
And d. Top Left: The semi-major axes evolve without perturbation 
for 200,000 years before planet d is suddenly flung from the system. 
Top Right: This configuration experiences low amplitude chaos in its 
eccentricity until it suddenly experiences a close approach at 200,000 
years. Bottom Left: The inclinations, however, evolve regularly until 
a close approach sends b into a 45° orbit. A second close approach 
kicks it back down, but d is ejected before it returns to the orbital 
plane. Bottom Right: This plot shows the typical behavior of A; it 
librates while motion is regular, but circulates after close approaches 
destabilize the system. Here, the A distribution represents 200,000 
years of libration and 25,000 years of circulation. 
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Fig. 31. — The long-term A distribution of the same 4 v And simula- 
tions as in Fig. 30. Top Left: This A distribution is typical of a regular 
system. Other Panels: The distributions are typical of chaotic config- 
urations. From these plots we conclude the secular resonance locking 
is important in maintaining stability, or regularness in v And. 

4.2. 47UMa 

The 47UMa system consists of a 1.03 ±0.03M Q (Gonzalez 
1998) star and 2 interacting companions: b and c, at 2.09AU 
and 3.73AU respectively. The initial eccentricities in this sys- 
tem are substantially lower than v And at 0.06 and 0.1 5 respec- 
tively. More recent measurements place e c much closer to 0. 
However, "e c =0.3 provides almost as good a fit to the radial ve- 
locity data as does e c =0.005" (Laughlin, Chambers, & Fischer 
2002). Should e c < 0.1 then, of all the systems examined in this 
paper, 47UMa most closely resembles our own. Planet b is the 
larger planet and hence determines T4 7( j Ma =840,000 years. 

Overall 80.3±4.7% of simulations were stable over t^wmo- 
This is less than a 2a difference from v And. In unstable con- 
figurations planet c, the less massive planet, was ejected every 
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Table 8 - Results of Long Term (10 8 yr) Simulations of v And 



Trial 


e c ,o 


ed,o 


Ro 


Ao (°) 


e 


Result" 


006 


0.302 


0.345 


5.351 


43.6 


2.1 x 10~ 9 


R 


007 


0.438 


0.371 


5.375 


7.8 


0.093 


C, e(2.6), E(b,2.7) 


054 


0.145 


0.201 


5.559 


56.6 


1 x 10~ 8 


C 


109 


0.392 


0.411 


5.28 


43.4 


0.32 


C, e(3.2), E(c,79.4) 


152 


0.253 


0.348 


5.378 


48.2 


2.1 x 10~ 8 


R 


164 


0.0938 


0.128 


5.422 


46.4 


9.5 x 10" 9 


R 


192 


0.252 


0.342 


5.448 


51.8 


3.9 x 10~ 5 


C, E(b,54.8) 


288 


0.348 


0.389 


5.227 


66.2 


1.2 x 10" B 


C 


499 


0.256 


0.349 


5.369 


1.1 


5.2 x 10~ 9 


R 


880 


0.242 


0.343 


5.138 


16.7 


1.6 x 10~ 8 


R 


989 


0.245 


0.305 


5.431 


26.0 


1.075 x 10- 4 


C, e(77.3) 



a R=Regular, C=Chaotic, e=Energy conservation failed (Time (10 6 yr)), P=Perturbed 
(Planet, Time (10 6 yr)), E=Ejected (Planet, (10 6 yr)) 
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time. This is similar to v And in which the most massive planet 
perturbed the smaller planets. The instability rate as a function 
of time is presented in Fig. 32. The rate for 47UMa is similar to 
the other systems in that most unstable configurations perturb 
a planet past stability on relatively short timescales, but with a 
tail to longer times. The rate does not reach zero, however, and 
suggests that more unstable configurations exist. 
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Fig. 32. — Instability rate in 47UMa. Most ejections occur at 10,000 
to 100,000 years. In this system, even unstable configurations gener- 
ally survive for at least 10,000 years. 



The 47UMa stability map is presented in Fig. 33. The over- 
all structure of stability in eccentricity space is qualitatively the 
same as in v And, with one major exception: stability is cor- 
related with ex (e c ), not e\ (e^). The errors for e c are substan- 
tial. The main difference between 47UMa and v And is that the 
more massive companion is the interior planet. This configura- 
tion makes it more difficult for the exterior planet to perturb the 
interior planet which is more tightly bound to the parent star. 

This stability map is in good agreement with other work done 
on this system. Using MEGNO, Gozdziewski (2002) deter- 
mined the maximum value for e c to be approximately 0.15. A 
stability analysis in Laughlin, Chambers, & Fischer (2002) also 
shows a similar dependence on e c . For nearly coplanar systems, 
such as those presented here, they determined the maximum 
value for e c to be less than 0.2. Although the exact maximum 
value for e c is different for all three studies, it is clear that the 
value of e c determines stability for 47UMa. 

When comparing Fig. 33 with Fig. 22, we see that the edge in 
the v And system is much cleaner than in 47UMa. One possible 
reason for this is the system's proximity to the 5:2 mean motion 
resonance. To examine this possibility, in Fig. 34 we plot stabil- 
ity as a function of R in this system. Although there appears to 
be some increase in stability beyond 5:2, and a decrease inside 
5:2, the errors are too large to confirm that this is a real effect. 
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Fig. 33.— The stability map for 47UMa. The relevant orbital ele- 
ments are ez, and e c . For this system the value of e c determines stabil- 
ity. The current best fit to this system, as of 1 Nov 2002, is marked 
by the asterisk (Laughlin, Chambers, & Fischer 2002). Note that al- 
though this system appears to lie far from the edge, the observational 
error for e c is ±0.1 15. 

In Figs. 35-37 we examine several possible orbital evolutions 
for 47UMa. Some sample simulations from the suite of 1000 
are listed in Table 9. First, in Fig. 35, is a stable regular con- 
figuration. Over 50% of all configurations for this system are 
regular. In this case the eccentricities, and inclinations show 
very small oscillations ( 10%) and A librates with an amplitude 
of 45°. 
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Table 9 - Selected Simulations of 47UMa 



Trial 


eb,o 


e c ,o 


Ro 


Comments" 


006 


0.021 


0.211 


2.384 


C, P(c,21.5), E(c,38.1) 


025 


0.043 


0.077 


2.503 


C 


032 


0.065 


0.105 


2.405 


R 



a R=Regular, C=Chaotic, P=Perturbed(Planet, Time (10 3 yr)), E=Ejected(Planet, 
Time (10 3 yr)) 
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Fig. 34. — The stability of 47UMa as a function of R. This sys- 
tem lies very close to the 5:2 mean motion resonance. There is some 
indication of more stability exterior to 5:2, and less interior, but the 
statistics are not robust enough to confirm this. The dashed line rep- 
resents the current best fit to the system. 
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Fig. 35. — The orbital evolution of 47UMa-032, a stable, regular 
configuration of 47UMa. The data are smoothed on 25,000 year inter- 
vals. Top Left: The semi-major axes show no hint of perturbations. Top 
Right: The eccentricities vary on a 4800 year timescale for the dura- 
tion of the simulation. Bottom Left: This system never deviates more 
than 5° from coplanarity. The inclinations oscillate on a 4400 year 
period. Bottom Right: In this configuration A is aligned, but librates 
with an amplitude of 45". 



As with all other systems, there are chaotic, but stable ex- 
amples of 47UMa. One such system, 47UMa-025, is exam- 
ined in Fig. 36. The semi-major axes of this configuration show 
no signs of perturbations, but the eccentricities and inclinations 
show clear signs of chaos, though no perturbation is very large. 
As with stable, chaotic configurations of other systems, circula- 
tion appears to be the dominant mode of A evolution. However, 
this system is actually librating about an anti-aligned orienta- 
tion, with only occasional circulation. The amplitude of this 
libration is quite variable, sometimes as high as 165°. 

In Fig. 37, we present an example of the ejection of planet 
c. This system evolves regularly for 10,000 years, then expe- 
riences 25,000 years of chaotic evolution, culminating in the 
ejection of the outer planet. For this configuration A circulates 
during the first 30,000 years, both regular and chaotic epochs, 
but with different angular velocities in each stage. During the 
final 10,000 years, however, when e c becomes very large, A 
becomes locked at 155°. 
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Fig. 36. — The orbital evolution of 47UMa-025, a stable, chaotic 
configuration of 47UMa. Top Left: The semi-major axes show no hint 
of perturbations. Top Right: The eccentricities vary by only 25% for 
the duration of the simulation, but show clear coupled, chaotic behav- 
ior. Bottom Left: This nearly coplanar simulation also shows obvious 
chaos in the evolution of the inclinations. Bottom Right: This A dis- 
tribution function is libration about anti-alignment with an amplitude 
which varies between 100° and 165°, but also with occasional circu- 
lation. 
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Fig. 37. — The orbital evolution of 47UMa-006, a chaotic, unstable 
configuration of 47UMa. The p marks the time that our perturbation 
criterion is met. Top Left: The semi-major axes evolve quiescently 
for 10,000 years, then e c begins its slow trek to upward of 200AU. 
The system is perturbed, however, until 21,500 years. Note that the 
y-axis is logarithmic in this example. Top Right: After 10,000 years, 
the system suddenly becomes chaotic, eventually pushing e c to unity 
in 40,000 years. Bottom Left: As with the eccentricities, the incli- 
nations evolve regularly for 10,000 years, but then become chaotic. 
Bottom Right: In this configuration A circulates, but eventually be- 
comes locked at 155° for the final 10,000 years. 
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Table 10 - Results of Long Term (10 9 yr) Simulations of 47UMa 



Trial 


eb,o 


e c ,o 


Ro 


Ao (°) 


e 


Result" 


238 


0.0587 


0.0809 


2.4006 


124.0 


4 x 10- 6 


R 


257 


0.0637 


0.0169 


2.3737 


57.7 


3.7 x 10- 6 


R 


307 


0.0594 


0.115 


2.4312 


122.6 


1.9 x 10" e 


R 


355 


0.0577 


0.0887 


2.4332 


64.3 


4.2 x 10" 6 


R 


507 


0.0607 


0.111 


2.4498 


40.8 


3.4 x 10" e 


R 


546 


0.063 


0.113 


2.3402 


6.1 


4.3 x lO" 6 


R 


624 


0.0595 


0.118 


2.5513 


164.5 


3.4 x 10- 6 


C 


633 


0.0583 


0.0900 


2.4872 


94.3 


4 x 10"° 


R 


792 


0.062131 


0.0897 


2.3957 


88.2 


2 x 10- 6 


R 



a R=Regular, C=Chaotic 
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Fig. 38. — Eccentricity evolution of 4 long term simulations of 
47UMa. Top Left: Evolution of 47UMa-624. Top Right: Evolution of 
47UMa-307. Bottom Left: Evolution of 47UMa-238. Bottom Right: 
Evolution of 47UMa-257. Most systems are regular for the duration 
of the simulation. However the top right is clearly chaotic, yet it still 
survives for 10 9 years. 

Long term simulations of 47UMa were integrated for 10 9 
years. Table 10 is a summary of all long term simulations for 
47UMa. Fig. 38 shows the eccentricity evolution of 4 systems. 
The top right panel of this figure is fascinating. The system is 
very chaotic for the first 350 million years, then enters a short 
(~20 million years) period of quiescence, only to return again 
to a similar chaotic state. The other configurations all appear 
regular. Fig. 39, plots the A evolution. In 47UMa no secular 
resonance locking occurs. However these plots confirm the re- 
sults of Laughlin, Chambers, & Fischer (2002). They show 
that for low values of e c (< 0.1) the system should librate in an 
aligned configuration, but above 0. 1 the system should be anti- 
aligned. The chaotic case, as expected, has a flat distribution 
function. 
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Fig. 39. — The distribution of A for 4 stable long term simulations of 
47UMa. The chaotic system (top left) shows a nearly flat distribution 
in A. This suggests that A is behaving chaotically as well. Two regular 
systems (top right, and bottom left) show libration about anti-parallel 
configurations, whereas the bottom right librates about A = 0. 

4.3. Gas Giants 

Perhaps the most interesting aspect of the new ESPS is their 
comparison to the SS. The procedure outlined in §2 permits a 
comparison of the ESPS and the SS. We must first, however, de- 
termine how to vary the initial conditions of the gas giants. As 
can be seen in Table 6, the "error" associated with each planet 
is arbitrary. We have given a spread in initial conditions that 
is roughly similar to the percent error as listed in the ESPS. 
For example the periods are allowed to vary by approximately 
10%, but the eccentricities have a standard deviation of 0.1 for 
all planets. This procedure will allow us to create a stability 
map for the SS, but will make a comparison of the probabilities 
of survival less meaningful. 

The outer SS consists of 4 gas giants located between 5.2 and 
40AU. The gas giants are on much more circular orbits than the 
ESPS planets (Saturn has the largest eccentricity at just over 
0.05). Because of the large semi-major axes, tss corresponds to 
3.32 x 10 6 years. Compared to known ESPS the gas giants are 
relatively low-mass planets. In fact Uranus and Neptune could 
not be detected via the Doppler method at the precision level 
currently achieved (see §4.4). 

Chaos in the SS is well documented (e.g. Sussman & Wis- 
dom 1988, Sussman & Wisdom 1992, Varadi et al. 1999, Lecar 
et al. 2001). In fact Varadi et al. (1999) show that the Jupiter- 
Saturn system lies very near chaotic regions. They vary the 
semi-major axes of these two planets slightly (less than 1 %) and 
find that this is enough to identify broad chaotic regions. Below 
we show that by enlarging this variation, the system moves into 
total instability; ejections are inevitable. 

For the gas giants 85.3%±4.3% of the trials were unper- 
turbed for 3.3 xlO 6 years. In Paper I, we integrated 32 gas 
giant configurations. Of these 81% survived. As in §4.1 we 
again recover the results of Paper I. In this system Jupiter was 
never ejected; it always removed the other planets from the sys- 
tem. Saturn was ejected in 14% of the simulations, Uranus, 
(the least massive) 62%, and Neptune, 24%. In the SS there- 
fore, ejection rate is tightly coupled to mass, as was observed 
in the other ESPS. It therefore seems that the SS behaves like 
other interacting systems. 

Fig. 40 is the instability rate. Most configurations survive 
for 10 5 years. Once again it appears the perturbation rate does 
not fall to zero, and we note that this means that we have not 
found all the unstable systems. The last bin in this plot contains 
only 20,000 years worth of data. It is therefore unclear if the 
ejection rate might still be rising with larger time. Should that 
be the case that our choice for t ss is too small, and suggests 
that we have not identified all unstable configurations. Fig. 41 
is the stability map for the gas giants. The gas giants show a 
plateau as in 47UMa and in v And, however the edge is much 
less dramatic. The actual values for our gas giants shows that 
our system lies quite far from the stability edge. We also note 
that the stability plateau shows many depressions, and the abyss 
contains many spires. This apparent difference between the SS 
and other interacting systems may result from not determining 
all unstable configurations. Perhaps longer integrations would 
sharpen the edge, and broaden the instability abyss. Dynami- 
cally speaking the major difference between the SS and other 
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interacting systems is that the SS has a much broader range of 
orbital times. Jupiter orbits 13 times more quickly that Nep- 
tune. Perhaps instability is more relevant on timescale based on 
the orbit of Neptune (see §6). However these features may also 
arise from the system's proximity to the 5:2 resonance, the so 
called "Great Inequality". We address this possibility in §4.4. 
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Fig. 40. — Rate of instability in the SS. Instability requires at least 
30,000 years to develop, and continues through tss- 
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Fig. 41 . — Stability map for the gas giants. In eccentricity space, the 
SS lies near a small depression. The edge in the gas giant system is 
not nearly as clean as in other interacting systems. This may because 
our choice of r is too low to find most unstable configurations. 



Several example simulations are shown in Figs. 42-47. The 
initial conditions and outcomes of these simulations are shown 
in Table 11. A regular, stable example is shown in Fig. 42. 
The eccentricities and inclinations are a linear superposition of 
eigenmodes. Although A begins nearly anti-aligned, it quickly 
moves into a large amplitude librational mode. The circulation, 
which is typically an indicator of chaos, appears to not influ- 
ence the eccentricity evolution. 

The best fit to ej, es and R is simulation SS-183. The orbital 
evolution of this system is shown in Fig. 43. Note, however, that 
Ao differs substantially from its standard value of 68.5°. This 
difference is responsible for the chaotic evolution of ej and e s . 
Curiously, though, the evolution of ejj and are regular. Un- 
like the eccentricities, all the inclinations evolve regularly. The 
nodes of Uranus and Neptune librate about anti-alignment, but 
with occasional circulation. Note that in e and ;' Uranus oscil- 
lates from 2 modes, whereas Neptune experiences 3 modes. 
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Table 11 - Selected Simulations of the Gas Giants 



Trial 


ej,o 


es,o 


Ro 


Comments" 


157 


0.155 


0.042 


2.40 


C, E(N,2.8) 


180 


0.155 


0.0.095 


2.45 


C, E(U,2.1) 


183 


0.056 


0.055 


2.48 


R/C 


278 


0.0944 


0.223 


2.45 


C, E(S,0.66) 


306 


0.124 


0.230 


2.46 


C 


402 


0.093 


0.203 


2.55 


R 



a R=Regular, C=Chaotic, E=Ejected(Planet, Time (10 6 yr)) 
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Fig. 42. — The orbital evolution of SS-402, a stable, regular configu- 
ration of the gas giants. The 4 planetary bodies in this system result in 
more complicated dynamics as the smaller 2 planets often experience 
several modes of oscillation. Top Left: The semi-major axes do not 
change for the duration of the simulation. Top Right: The eccentric- 
ities are combinations of sinusoids. Jupiter and Saturn both oscillate 
on a 70,000 year period. Uranus and Neptune vary on a 1.5 million 
year timescale. Their motion however also has 4 million year mode, 
due to to beating with the oscillations from Jupiter and Saturn. Bottom 
Left: All the inclinations evolve as a function of two modes. Jupiter 
and Saturn oscillate on a 40,000 year cycle. Uranus and Neptune os- 
cillate on a 500,000 year period. As with eccentricity we observe 
beating between these 2 frequencies in all planets. Bottom Right: A 
switches between libration with a 80° amplitude, and circulation. 
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Fig. 44. — The orbital evolution of SS-306, a stable, chaotic con- 
figuration of the gas giants. Top Left: The semi-major axes begin 
migrating immediately, but a change greater than 20% does not oc- 
cur for the duration of the simulation. Top Right: All eccentricities 
undergo chaotic oscillations, but the amplitudes are small. No eccen- 
tricity ever reaches 0.35. Bottom Left: The inclinations also evolve 
chaotically with low level oscillations. The system remains close to 
coplanarity, as is never exceeds 12°. Bottom Right: This A distribu- 
tion is clearly chaotic as the 2 peaks and circulation demonstrate. 



In Fig. 45 we present an example of the ejection of Saturn. 
Although the semi-major axes show obvious signs of perturba- 
tions, it is only during the last 50,000 years that any planet's 
semi-major axis changes by more than 10%. The eccentricities 
and inclinations do not show large fluctuations until Saturn is 
quickly ejected after 65,000 years. 



Fig. 43. — The orbital evolution of SS-183, a stable configuration of 
the gas giants in which some elements evolve regularly, others chaot- 
ically. Top Left: The semi-major axes do not change for the duration 
of the simulation. Top Right: The eccentricities of Jupiter and Sat- 
urn evolve chaotically, but Neptune and Uranus appear to be regular. 
Bottom Left: All inclinations evolve regularly, although the number 
of modes is different, Jupiter, Saturn, and Neptune have 2 modes, but 
Uranus has 3. Bottom Right: A shows the typical distribution func- 
tion of libration about anti-alignment and circulation. This alternating 
results in the chaotic evolution of ej and es through e - 05 coupling. 



In Fig. 44 a fully chaotic, yet stable configuration is shown. 
The semi-major axes show obvious signs of encounters, but 
never change by more than 20%. The eccentricities are very 
chaotic, but rarely reach 0.3. The inclinations, too, are ex- 
tremely chaotic, but never surpass 12°. The double peaked A 
distribution function is another clear indicator of chaos. 
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Fig. 45. — The orbital evolution of SS-278, the ejection of Saturn. 
Top Left: The semi-major axes deviate from the initial values by less 
than 10% until the final 50,000 years of integration. Even then as only 
increases by 15%. Top Right: All eccentricities are chaotic, but the 
values remain close to their initial conditions until Saturn is thrown 
from the system. Bottom Left: The inclinations appear to evolve regu- 
larly for 35,000 years until Jupiter and Saturn become chaotic. Uranus 
becomes chaotic after 60,000 years, but Neptune remains regular for 
the duration of the simulation. Bottom Right: The A distribution for 
this configuration is quite unusual and also indicates the system is 
chaotic. 

The ejection of Uranus is shown in Fig. 46. The semi-major 
axes of Jupiter and Saturn remain nearly constant throughout 
the simulation, but Neptune and Uranus are clearly interact- 
ing. The eccentricities are chaotic, but remain near the starting 
values, except for Uranus, which steadily increases until it is 
ejected after 2 x 10 6 years. The inclinations are also chaotic. 
Jupiter and Saturn appear to have 1 normal mode, and 1 chaotic 
mode, while Neptune and Uranus are fully chaotic. The lon- 
gitudes of periastron tend to remain near alignment, but the 3 
peaks clearly bely the chaos in the system. 
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Fig. 46. — The orbital evolution of SS-180, the ejection of Uranus. 
Top Left: The semi-major axis of Jupiter does not change, and Saturn 
changes by less than 0.1 AU at the very end of the simulation. Uranus 
and Neptune are clearly interacting, but the fluctuations of eu and E N 
are small until 1.75 x 10 6 years. Top Right: While the eccentrici- 
ties of Jupiter, Saturn, and Neptune remain low, Uranus's eccentricity 
gradually grows until it is ejected after 2 x 10 6 years. Bottom Left: 
The inclinations all also show chaos, but Jupiter and Saturn appear to 
have a regular 3000 year mode superimposed on small chaotic fluc- 
tuations. Bottom Right: The A distribution for this configuration is 
quite unusual and also indicates the system is chaotic. Note, however, 
that the system is experiencing a generic form of libration, as A never 
exceeds 60°. 

Finally we plot the evolution of a configuration that ejects 
Neptune in Fig. 47: SS-157. After 750,000 years Uranus be- 
comes the most distant planet. After Uranus and Neptune have 
a close approach which flings Neptune from the system, Uranus 
drops to a 15AU orbit. It is likely that it would then encounter 
Saturn, and another planet would be removed from the system. 
Despite the large variations in semi-major axes, the eccentric- 
ities remain relatively calm. Jupiter and Saturn in particular 
appear nearly regular. In inclination, they in fact are regular. 



For all the unstable cases shown here, this configuration's A 
shows the least chaos. Although there is some flipping between 
a librational and a circulation mode. A also demonstrates that 
the chaos between Neptune and Uranus does not affect the evo- 
lution of Jupiter and Saturn's orbits. 

These last three figures reveal why the stability map for the 
SS is different than for other interacting systems. Sometimes, 
even when ej and e s are low, the smaller bodies in the system 
can interact with each other violently. Although the eccentricity 
of Jupiter and Saturn are generally the most important factor in 
system stability. The larger the values of any eccentricity, the 
more likely is instability. 

We ran no long term simulations on the SS. The orbital el- 
ements for the SS are well determined, and many long term 
simulations have already been performed on this system (see 
Duncan & Quinn 1993, Lecar et al. 2001). 
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Fig. 47. — The orbital evolution of SS-157, the ejection of Nep- 
tune. These data are smoothed over a 40,000 year interval. Top Left: 
After 750,000 years Neptune and Uranus cross orbits, and Uranus be- 
comes the most distant planet. After ejecting Neptune Uranus drops 
to a 15AU orbit which would most likely encounter Saturn relatively 
quickly and result in another ejection. Jupiter and Saturn remain on 
regular orbits. Top Right: All eccentricities remain low for the dura- 
tion of the simulation until Neptune is removed after 2.8 x 10 6 years. 
This is true even while Uranus and Neptune are in 1:1 resonance at 
750,000 years and 2.2 million years. Jupiter and Saturn experience 
sinusoidal oscillations throughout the simulation. Bottom Left: The 
inclinations of Jupiter and Saturn are regular but Uranus and Neptune 
are clearly chaotic. Note that the encounter that throws Neptune from 
the system sends Uranus into a highly inclined orbit. Bottom Right: 
For this configuration A is generally librating in an aligned state with 
an amplitude of 70°, but with some circulation. 

4.4. Jupiter and Saturn 

As mentioned above Uranus and Neptune do not provide 
enough reflex velocity, K, motion in the star to be observable 
by current technology (K ~ 3m/s). Should any planet of Uranus 
or Neptune mass exist in the observed ESPS they would not be 
detected. Therefore we followed the same procedure with just 
Jupiter and Saturn. This suite of simulations can also provide 
clues as to how other ESPS will behave if they have additional, 
distant companions. 
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Not surprisingly this 3-body system is more stable than the 
5-body system, as 96.3 ±2.4% of the trials remained stable. In 
this system Jupiter was ejected 6.9% of the time, and Saturn 
93.1%. It seems as though instability is being passed through 
Saturn and into the smaller planets. We can therefore apply this 
result to the other ESPS. This Jupiter-Saturn system is analo- 
gous to the 47UMa system. The mass ratio of the planets are 
about the same, as is R. The only substantial difference is that 
the masses are higher and the orbits closer in 47UMa. However, 
our simulations of 47UMa used e c =0. 1 . This is about twice as 
high as e Saturn . This decrease in eccentricity is clearly impor- 
tant as the Jupiter-Saturn system is substantially more stable 
than 47UMa. This again demonstrates that the eccentricities of 
the planets are key in determining stability. 




Fig. 48. — Stability map for the Jupiter-Saturn system. This map 
shows some similarities to that for the gas giants (see Fig. 44). The 
line demarcating the plateau follows approximately the same diagonal 
line. As well we still see a few depressions in the plateau. 

Fig. 48 is the stability map for the Jupiter-Saturn system. 
This plot is similar to Fig. 41. There is a boundary at approx- 
imately the same location in eccentricity space, however the 
drop is not so sharp, nor as deep as in the gas giant system. 
Further an additional plateau rises at larger eccentricities. This 
last phenomenon is not observed in any other system in this 
paper. 

As mentioned in §4.3, stability might be correlated with the 
5:2 resonance. Fig. 49 shows stability as a function of R for the 
SS. There is a hint that as the system moves out of this third 
order resonance, stability increases, but the errors on these dis- 
tant configurations are too large to confirm this possibility. In 
47UMa R = 2.38, while in the SS, R = 2.48. Although there are 
no statistically meaningful points in Figs. 34 and 49, the same 

6 http://obswww.unige.ch/~udry/planet/hd83443_syst.html 
7 http://obswww.unige.ch/~udry/planet/hd74156.html 



trend appears in both. Namely that interior to 5:2 instability is 
more prevalent. As has been shown throughout this paper, the 
eccentricities determine the overall stability, and the statistics 
are too poor to claim any trend with R in either system. 

Comparing the gas giants with Jupiter-Saturn provides us 
with an excellent opportunity to explore completeness. As men- 
tioned above, the Jupiter-Saturn system is similar to 47UMa, 
yet the stability maps are quite different. The Jupiter-Saturn 
system is the only system examined with no instability abyss, 
and it is the only system we know to be incomplete. 




Fig. 49. — Stability as a function of R for the Jupiter-Saturn sys- 
tem. As in 47UMa there is a hint that instability increases inside the 
5:2 resonance; no data point lies more than la from the mean rate 
of stability (96%). It therefore apears that this resonance has mini- 
mal impact on the system. The dashed vertical line represents the true 
value of R in the SS. 

5. SEPARATED SYSTEMS 

By the end of 2002 three separated systems had been an- 
nounced: HD83443 6 , HD168443 (Marcy et al. 2001b), and 
HD74156 7 . The values and errors for these systems are re- 
produced in Table 12. HD83443 consists of 2 Saturn mass 
planets in very tight orbits. HD168443 consists of 2 very 
large companions im\ > \!M Jup , ni2 > 1 .5M Jup ). In fact planet 
c should be considered a brown dwarf, and if the system is 
more inclined than 35 degrees planet b would also be a brown 
dwarf. For this system R = 30.5. HD74156 contains 2 bodies of 
slightly more than a Jupiter mass, with R = 44.6. We only ex- 
amined HD168443 and HD83443. Evidence is mounting that 
HD83443 is not a multiple system (Butler et al. 2002), therefore 
we stopped the simulations on this system after 847 trials had 
been completed. For HD168443 and HD83443, all simulations 
survived to r. HD74156 has a larger/?, and smaller masses than 
HD 168443. It therefore seems highly doubtful that any simula- 
tion of HD74156 would produce an unstable configuration. 

Although all simulations were stable, the dynamics of 
HD 168443 are still interesting. The eccentricities and inclina- 
tions of this system show a weak planet-planet interaction. Al- 
though no evidence of chaos is evident, the planets apparently 
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Table 12 - Initial Conditions for Separated Systems 



System 


Planet 


Mass(Mj) 


Period(d) 


Eccentricity 


Long, of Peri.(Deg) 


Time of Peri. (JD) 


HD168443 


b 


7.73 


58.1 ± 0.006 


0.53 ±0.003 


172.9 ±0.4 


2450047.58 ±0.2 




c 


17.15 


1770 ± 25 


0.20 ±0.01 


62.9 ±3.2 


2450250.6 ± 18 


HD74156 


b 


1.56 


51.619 ±0.053 


0.649 ±0.022 


183.7 ±3.3 


2451981.4 ±0.57 




c 


7.3 


2300(Fixed) 


0.395 ±0.074 


240 ± 12 


2450819 ±75 


HD83443 


b 


0.34 


2.9853 ±0.0009 


0.079 ±0.033 


300.25 ± 17.05 


2451386.5 ±0.14 




c 


0.16 


29.83 ±0.18 


0.42 ±0.06 


337.42 ± 10.42 


2451569.59 ±0.73 
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are close enough that they feel each other. This system may 
be fully stable, but it appears to lie close to the boundary be- 
tween interacting and separated systems. We hypothesize that 
this proximity to interacting systems is due to the large plane- 
tary masses in this system. 

6. SUMMARY 

In this paper we have described the dynamics of three differ- 
ent morphological classifications of planetary systems. We find 
that the systems in each of the classifications have similar sta- 
ble regions. In resonant systems very small stability zones exist 
in phase space, and stability is tightly coupled with R and, to a 
lesser degree, e\, where e\ is the eccentricity of the most mas- 
sive companion. In interacting systems, the zones are larger, 
but are correlated with e\ and ei, where is the eccentricity of 
the second most massive companion. In these systems we see a 
correlation with eccentricity and the location of the most mas- 
sive planet. Large interior planets are almost impossible to eject 
(47UMa, the gas giants), whereas large exterior planets can be 
ejected sometimes (v And). Separated systems are completely 
stable as observed. 

Table 13 summarizes the results of this paper. In this ta- 
ble /stable is the fraction of configurations that were stable, and 
fjj = 1,2,3, ... is the fraction of unstable systems which per- 
turbed/ejected that planet. Note that the subscripts correspond 
to mass, not semi-major axis (1 being the most massive com- 
panion). We therefore strengthen the theory suggested in Paper 
I: all interacting planetary systems lie near the edge of stability. 

There are some obvious similarities among the systems. One 
particularly intriguing result is the similarity in f stable between 
systems in the same classification. Resonant systems have sur- 
vival probabilities less than 20%, whereas interacting systems 
lie close to 80% and the separated systems are completely sta- 
ble. Of the interacting systems, 47UMa stands out as being far 
from the stability edge. 

The choice of r = 2.8 x 10 5 Pi appears to identify most un- 
stable configurations. For both resonant systems r < 10 6 yrs, 
but since the simulations were integrated to 10 6 yrs we may es- 
timate the usefulness of this arbitrary value. For both GJ876 
and HD82943 approximately 2-3% of configurations ejected a 
planet after r. For the v And system, 1.5% of unstable sys- 
tems were ejected in the last bin of Fig. 20. For 47UMa the 
rate was over 10%, and for the SS the rate was 4%. However 
as was noted in §4.3, this low rate for the SS may be the result 
of poor sampling in the last bin. Although all these systems 
reached a maximum ejection rate before r, the nonzero rate at 
t demonstrates that our choice for r was slightly too short. In 
GJ876 some ejections occurred right up to 10 6 yrs, but by 10 4 
yrs (0.25T G yg 7 6) over 90% of unstable cases had been identified. 
The situation is nearly the same for HD82943; 90% of unstable 
cases were identified by 30,000 years (0.1 77/082943), but ejec- 
tions continued for 10 6 years. Given these statistics, a better 
choice for t would be r = 10 6 P„ uter . However, it is important to 
note that instability can arise after this, as is shown in Figs. 9, 
18, and 30. The simulations presented here clearly demonstrate 
the unpredictable behavior of chaotic systems; no choice of t 
would identify all unstable configurations. One should there- 
fore note that all the global results presented here are upper 
limits. The probability of survival and extent of stable phase 
space are smaller than what is shown here. 

Long term simulations (> 10 8 orbits) show that all systems 
have regular configurations on this timescale. However only 
one simulated configuration of GJ876 showed this behavior. 



Some configurations may show a large degree of chaos for up 
to 10 9 years (see Fig. 38, top left), eject a planet after an arbi- 
trarily long period of time (Fig. 30, bottom right), in addition 
to quiescent, regular evolution. Regular orbits tend to librate 
about A=0, but this is not necessary for stability (see Figs. 10, 
19, 31, and 39). 

7. DISCUSSION 

Ideally this research provides insights into planet formation. 
In particular, the current distribution of orbits may give clues to 
the formation scenario. Two features of ESPS are particularly 
interesting: the apsidal alignments and the large eccentricities. 
As e and 03 are coupled, these two phenomena are likely the 
result of the same mechanism. There are two generic ways to 
pump up eccentricities: adiabatically, or impulsively, with re- 
spect to the secular timescale (;> 10 5 yr). Our variation of or- 
bital elements provides a unique view into the effects of these 
mechanisms on the dynamics and stability of actual planetary 
systems. Several groups have examined this problem, and in 
this section we interpret our results in the context of theirs. 

Of adiabatic scenarios a remnant planetary disk is the most 
likely candidate (Chiang & Murray 2002; Goldreich & Sari 
2003). For at least the v And system, a remnant disk exter- 
nal to planet d can provide a mechanism to pump up e c and 
to their current observed values (Chiang & Murray 2002). This 
method also predicts libration of A about 0, which is observed 
in this work. Conversely an impulsive force may also drive ec- 
centricities to values significantly higher than zero (Malhotra 
2002). The impulsive scenario also perturbs A. In adiabatic 
schemes the libration amplitude is small, whereas in impulsive 
cases it can be quite large (> 45"). Throughout this paper we 
have shown configurations with libration amplitudes larger than 
45°, i.e. Figs. 5 and 39, and smaller than 45°, i.e. Figs. 10, 31. 
This work therefore finds examples of systems which may re- 
sult from either mechanism. 

This impulsive scenario is difficult to reconcile with reso- 
nant systems. These systems most likely form as a result of 
resonance capture during the orbital migration epoch of planet 
formation (Snellgrove, Papaloizou, & Nelson 2001), which as- 
sumes adiabatic migration. This phenomenon seems qualita- 
tively similar to the external disk model of Chiang & Murray. 
Their model, based on torques produced by Lindblad and coro- 
tation resonances, is very similar to planetary migration. How- 
ever current orbital migration theory predicts that a Jupiter mass 
planet at 5 AU in a plausible minimum mass solar nebula should 
migrate on a timescale of order 2500 years (Tanaka, Takeuchi, 
& Ward 2002) to 5000 years (Lufkin et al. 2003), which is a 
factor of 5 to 10 times shorter than the typical secular timescale 
for planetary systems. This suggests the planetary disk model 
might actually be impulsive, but only marginally so. However 
Lufkin et al. also point out that the migration might be very 
impulsive in heavier disks. Understanding the rates of migra- 
tion will be a major step toward resolving this issue of high 
eccentricities and apsidal alignment. All we can say now is that 
we are too limited in the number of resonant and interacting 
systems to determine if their eccentricities result from similar 
processes. 

The results presented here, coupled with Malhotra (2002) 
support the theory that the eccentricities of planets in interact- 
ing systems result from planet-planet scatterings. This possi- 
bility has been investigated substantially (Rasio & Ford 1996; 
Ford, Havlickova, & Rasio 2001; Weidenschilling & Marzari 
2002; Malhotra 2002). The proximity of these systems to the 
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Table 13 - The Stability of Planetary Systems 



System 


R 


T (10V) 


Ngood 


f stable 


h 


h 


h 


h 


GJ876 


2.03 


0.47 


950 


0.056 


0.04 


0.96 






HD82943 


2.01 


3.4 


955 


0.188 


0.09 


0.91 






47UMa 


2.38 


8.4 


997 


0.803 


0.00 


1.00 






Gas Giants 


2.48 


33.2 


996 


0.857 


0.00 


0.14 


0.24 


0.62 


Jup-Sat 


2.48 


33.2 


992 


0.963 


0.07 


0.93 






v And 


5.4 


10 


996 


0.861 


0.07 


0.54 


0.39 




HD168443 


30.5 


13.6 


1000 


1.0 


n/a 


n/a 






HD83443 


9.99 


0.023 


847 


1.0 


n/a 


n/a 
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edge of stability might imply that planet formation is an effi- 
cient process. Perhaps too efficient. As planets form, they are 
constantly perturbing each other with ever greater force. It is 
well known that ejections are common during planet formation. 
In fact some research predicts that the ejection of a fifth terres- 
trial planet may be needed to explain the period of heavy bom- 
bardment in the SS (Chambers, Lissauer, & Morbidelli 2001). 
So it is not too surprising that we find systems near instability 
because they form in an unstable state and eject massive bod- 
ies until they arrive in resonance, or reach the stability plateau. 
Some work has shown that if planet formation is very efficient 
(i.e. initially 10 Jupiter mass planets) then the subsequent scat- 
tering and ejections can produce distributions of a and e that are 
similar to those observed (Adams & Laughlin, 2003). Clearly 
this scenario is appealing, and will be verified in the next sev- 
eral years as simulations become more sophisticated and more 
multiple planet systems are detected. 

Beyond the origin of large eccentricities and apsidal align- 
ment, we find some inconsistencies in the theory of the origin 
of the very short period (P < lOdays) planets. As mentioned in 
§2 and §4. 1, the effects of tidal circularization were not included 
in this model. For v And the timescale for circularization is 80 
million years (Trilling 2000), but we see that the eccentricity 
of v And b can oscillate on 10 5 year timescales with an am- 
plitude of 0.3. At this point is is unclear if the tidal damping 
will always overwhelm the perturbations of other companions. 
It could be that we have detected v And b at a point in time 
in which its orbit is nearly circular. Perhaps we will discover 
planetary systems in which a close planet is being perturbed by 
external companions and the tidal circularization cannot com- 
pensate. However it seems more likely that the circularization 
is a stronger effect as we have yet to detect any planet inside 
the circularization radius on an eccentric orbit. Future numeri- 
cal work should resolve this issue. 

We observe that, in general, there is good agreement between 
our results and those performed with MEGNO. The sizes of sta- 
ble regions appear to be overestimated in those papers, but that 
is most likely due their choice for r, usually 10 4 years. The 
shortcomings of MEGNO are most clearly demonstrated in the 
long term integrations of 47UMa. The top left of Fig. 38 shows 
a stable, but chaotic system which persists for 10 9 years. The 
uniqueness of systems such as this is unknown, but understand- 
ing the dynamics of chaotic, yet long-lasting, systems could 
yield new insights into planetary dynamics. Our own SS is an- 
other example of a system that displays weak chaos, yet can 
survive for very long periods of time (Laskar 1994). 

This research is the first to examine the origin of high eccen- 
tricities and apsidal alignment with known systems. Most other 
work in this field is purely hypothetical. That type of research 
has the benefit of being unconstrained by statistics; they may 
integrate as many systems as they wish, with arbitrary initial 
conditions. This work, conversely, is the first attempt to co- 
herently and consistently compare known systems in order to 
understand their dynamics and origins. At this point, with so 
few known systems, the two methods are complimentary, but 
as we discover more ESPS, the method described in this paper 
will become more valuable as it uses true ESPS as a starting 
point. 

8. CONCLUSIONS 

We have shown that this type of experiment can indeed con- 
strain the observed orbital elements of planetary systems. Fur- 
ther we see that in almost all interacting and resonant systems 



the current best fits to the system place them near the bound- 
ary between stability and instability. The fact that no system is 
completely unstable implies that the observations of these sys- 
tems are reliable, and that the errors in the system are probably 
conservative. That is, all systematics have been removed, and 
statistical fluctuations are being overestimated. 

Note that our estimates of the instability of these systems is 
in some sense an underestimate because of the possible pres- 
ence of yet undetected lower mass companions. For example, 
it may turn out that 47UMa may have an undetected planet that 
would put it closer to the edge. On the other hand, the very 
existence of these systems shows that they are not unstable. As 
unsettling as it may be, it seems a large fraction of planetary 
systems, including our own, lie dangerously close to instabil- 
ity. As more and different types of systems are detected we will 
discover if all planetary systems are on the edge. 

This method has shown that, dynamically, the SS is not a 
unique system. In fact, it lies in the middle stability category. 
Some systems lie nearer instability, others further. As the radial 
velocity searches continue, and astrometric searches begin, a SS 
analogue (circular orbits, large semi-major axes) will undoubt- 
edly be discovered and we will finally be able to determine how 
the SS fits in with other planetary systems. But this experiment 
has shown that, with regards to its (close) proximity to unstable 
regions, the SS is a typical planetary system. 

Recently more systems were announced; HD38529 (Fischer 
et al. 2003), a separated system, HD 12661 (Fischer et al. 2003), 
an interacting system, and 55Cnc a 3 planet system with inte- 
rior planets in 3:1 resonance and a distant companion (Marcy 
et al. 2002). The planets in HD38529 have masses less than 
HD168443, and comparable values for R, therefore it seems 
likely that they are fully stable. 55Cnc, however, might demon- 
strate different dynamics and edges as it is in a different mean 
motion resonance. Future planetary systems will most likely 
fall into the categories outlined in this paper. The results pre- 
sented here suggest that f sla bie for 55Cnc would lie between res- 
onant and coupled systems. HD12661 is very similar to v And, 
so we expect this system to show similar edges, probabilities, 
and dynamics. 

Future work will address many of the issues brought up in 
§7. If planet formation is an efficient phenomenon, then we 
might suspect that additional companions lie in separated sys- 
tems. Further we should be able to determine that the eccentric- 
ity distribution of ESPS result from a late scattering event. We 
also need to determine the origin of the edges presented here. 
A mathematical relationship probably exists which would make 
the classification of planetary systems trivial. The categories as 
defined here may be the result of small number statistics. Two 
systems, HD169830 and HD37124, in which R w 10 have been 
announced (Mayor et al. 2003, Butler et al. 2003). These sys- 
tem may reveal the boundary between interacting and separated 
systems. An analysis of these two systems and 55Cnc will help 
sharpen our classification of planetary systems. 

Future work, both observational and theoretical, must ad- 
dress these issues. These systems as they are observed now 
reflect their histories, and hence provide us with the best path 
to unlocking the secrets of planet formation. As more and 
more observations of these planetary systems, additional plan- 
etary systems, and (hopefully) protoplanets are made, numeri- 
cal studies such as this, and those cited here, should provide a 
deeper understanding of planet formation and dynamics. 

9. ACKNOWLEDGMENTS 



Barnes & Quinn 



37 



The authors wish to thank Chance Reschke for his help in 
completing these simulations. This manuscript was greatly 
clarified after an edit by Derek C. Richardson. We would also 
like to thank Renu Malhotra, and Nader Haghighipour for use- 



ful discussions and suggestions. This work was funded by 
grants from NASA, NAI, and the NSF, and was simulated on 
computers donated by the University of Washington Student 
Technology Fund. 



REFERENCES 



Adams, F.C., & Laughlin, G. 2003, Icarus, 163, 290. 
Barnes, R. & Quinn, T. 2001, ApJ, 550, 884 
Benedict, G.F. et al. 2002, ApJ, 581, L115 
Butler et al. 1999, ApJ, 526, 916 
Butler et al. 2002, ApJ, 578, 565 
Butler et al. 2003, ApJ, 582, 455 

Chambers, J.E., Lissauer, J.J., Morbidelli, A. 2001, DPS, #33.2309 

Chiang, E.I., Tabachnik, S., Tremaine, S. 2001, AJ, 122, 1607 

Chiang, E.I., Murray, N. 2002, ApJ, 576, 473 

Cincotta, P.M., & Simo, C. 2000, A&AS, 147, 205 

Duncan, M.J., Quinn, T., 1993 ARA&A, 31, 265 

Fischer et al. 2002, ApJ, 564, 1028 

Fischer et al. 2003, ApJ, 586, 1394 

Ford, E.B., Havlickova, M., Rasio, F.A. 2001, Icarus, 150, 303 

Gonzalez, G. 1998, A& A, 334, 221 

Gonzalez, G., & Laws, C. 2000, AJ, 119, 390 

Goldreich, S., & Sari, R. 2003, ApJ, 585, 1024 

Gozdziewski, K. 2002, A&A 393, 997 

Gozdziewski, K. & Maciejewski, A.J. 2001, ApJ, 563, L81 

Konacki, M. & Maciejewski, A.J. 1999, ApJ, 518, 442 

Laughlin, G., Adams, F.C. 1999, ApJ, 526, 881 

Laskar, J. 1994, A&A, 287, L9 

Laughlin, G., Chambers, J., & Fischer, D. 2002, ApJ, 579, 455 
Lecar, M., Franklin, FA., Holman, M.J., & Murray, N.W. 2001, ARA&A, 39, 
581 



Levison, H.F, & Duncan, M.J. 1994, Icarus, 108, 18 

Lissauer, J.J. 1993, ARA&A, 31, 129 

Lufkin, G. et al. 2003, MNRAS, in press 

Malhotra, R. 2002, ApJ, 575, L33 

Marcy, G.W. et al. 2001a, ApJ, 556, 296 

Marcy, G.W. et al. 2001b, ApJ, 555, 418 

Marcy, G.W. et al. 2002 ApJ 581 1375 

Marzari, F., Weidenschilling, S.J. 2002, Icarus, 156, 570 

Mayor, M. et al. 2003, A&A accepted, astro-ph/0310316. 

Michtchenko, T., & Ferraz-Mello, S. 2001, AJ, 122, 474 

Rasio, FA., Tout, C.A., Lubow, S.H., Livio, M. 1996, ApJ, 470, 1187 

Rasio, F, & Ford, E. 1996, Science, 274, 954 

Rivera, E.J., & Lissauer, J.J. 2000, ApJ, 530, 454 

Rivera, E.J., & Lissauer, J.J. 2001a, ApJ, 558, 392 

Rivera, E.J., & Lissauer, J.J. 2001b, ApJ, 554, 1141 

Robutel, P., & Laskar, J. 2001, Icarus, 152, 4 

Santos, N.C., Israelian, G., Mayor, M. 2000, A&A, 363, 228 

Snellgrove, M.D., Papaloizou, J.C.B., Nelson, R.P. 2001, A&A, 374, 1092 

Stepinski, T.F., Malhotra, R., Black, D.C. 2000, ApJ, 545, 1044 

Sussman, G.J., & Wisdom, J. 1988, Science, 241, 433 

Sussman, G.J., & Wisdom, J. 1992, Science, 257, 56 

Tanaka, H., Takeuki, T., & Ward, W.R. 2002, ApJ, 565, 1257 

Trilling, D.E. 2000, ApJ, 537, L61. 

Varadi, F, Ghil, M., & Kaula, W.M. 1999, Icarus, 139, 286 



